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ABSTRACT 


A theoretical study is made of the visible and UV line radiation of 
He I atoms and He II ions from a plane -parallel model flare layer. Codes have 
been developed for the solution of the statistically steady state equations 
for a 50 level He I - II - III model, and the line and continuum transport 
equations. These codes are described and documented in the report along -with 
sample solutions. Optical depths and some line intensities are presented for 


a 1000 km thick layer. Solutions of the steady state equations are presented 

4 4 0 

for electron temperatures 10 - 5 x 10 K and electron densities 
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HELIUM EMISSION FROM ACTIVE SOLAR REGIONS 
I. INTRODUCTION 

The purpose of this program has been to develop codes for the 
simultaneous calculation of He I and II resonance line and He I D^ line 
intensities from model flare regions. These lines were chosen because 
of the spectral ranges of the Sky lab high resolution spectrographs and 
because of the planned program to obtain filtergrams on a patrol 
basis at the Lockheed Itye Canyon Observatory during the ATM mission. 

The NRL spectroheliograms incorporate simultaneous measurements of the 
He I and II resonance lines. 

A plane-parallel layer irradiated on one side by the photospheric 
radiation field was chosen as the geometric model. A statistically 
steady state and uniform electron temperature and density with position 
were assumed. The energy level model consists of all terms through prin- 
cipal quantum number 4. Our study has been confined to conditions we 

believe characteristic of flare regions, namely electron temperatures be- 

4 4 cl 10 14 -3 

tween 10 and 5 x 10 & and electron densities between 10 and 10 cm . 

An extensive compilation of electron impact excitation rates has been 

made as part of this study. The results were published in Solar Physics 

(Benson and Kulander, 1972). 

The statistically steady state level populations of model He I 
atoms have been calculated by a number of investigators for temperatures 
and densities characteristic of the outer solar atmosphere. Almost none 
of these authors has considered a sufficiently detailed energy level 
structure in the model atom to accurately obtain the line emission. 
Jefferies (1955) e.g. treats the 2s and 2p levels as a single level. 

De Jager and de Groot (1957) consider the 2s and 2p terms separately 
but the term structure of higher levels is ignored. This higher term 
structure is also ignored by Athay and Johnson (i960). They also neglect 
the effect of the He II ion processes by using other values for the 
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He i/He XI equilibrium. Athay and Johnson arrive at the conclusion that 

in the temperature range 40,000 - 50>000°K, the Do line will appear in 

12 ~ l 
emission for n g 10 almost independently of T . 

Zirin (1956) assumes in his calculations that transitions between 
terms of a given level are of negligible importance in determining the 
occupation numbers. This is known to be a poor assumption. Shklovsky 
and Kononovitch (1958) have calculated the D, line intensity but have made 
a number of unrealistic physical assumptions. More recently Hearn 
(1969) has calculated the occupation numbers of a 41 level He I atom and 
one level He II ion but he only presents results for the resonance line 
intensities. We shall demonstrate that more levels are required in He II 
to obtain correct line intensities. 

Jefferies (1957) has calculated the line intensity from a layer 
assumed to be optically thick in the D^ line. The transport equation was 
solved assuming incoherent scattering with no photospheric radiation in 
the line. Jefferies* results are very qualitative since it is known that 
the D^ line is probably not optically thick. 

To obtain accurate line intensities, simultaneous solution of the 
line and continuum radiative transfer equations and the steady state 
populations is required. To accomplish this, we have developed three 
separate codes. The first code (Code l) solves the statistical equili- 
brium equations for a 30 level He I-II-III system given the appropriate 
rates. The basic equations and sample solutions are given in Section II. 
The code is described in Section V and a listing is given in Appendix A. 

The second code (Code 2) represents a numerical solution of the 
line transport equations for a finite layer. The solution is of the 
integral form of the transport equation by expansion of the source func- 
tion in terms of a finite sum. The mathematical method used Is summarized 
by Avrett and Loeser (1969). Complete frequency redistribution and a 
Gaussian absorption profile are assumed. The basic equations and sample 
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solution are given in Section III. The code is described in Section V 
and a listing is given in Appendix B. 

The third code (Code 3) solves the continuum transport equation 
by expansion of the source function in a very similar manner to the line 
transport equation. The equations are given in Section XV. The code 
is described in Section V and a listing is given in Appendix C. 
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II. THE STEADY STATE EQUATIONS 

A. Energy Level Model 

The 50 assumed energy levels for the system of ions He I - He III are 
given in Table II -1. There are 19 levels for He I, 10 for He II and 1 for 
He III. These levels are shown in Figures II -1 and II-2. The levels are 
numbered 1 - 30 in order of increasing energy. Levels 1, 20 and 30 are the 
ground states of He I, II and III, respectively. This model was chosen to 
provide accurate solution for the first two resonance lines in He I and II 
and the D3 and IO830 lines of He I. Many other lines are included but were 
not the primary lines under consideration. The allowed transitions included 
in the model are listed in Table II. 2 together with f numbers and wavelengths. 
The inclusion of the 4 S,P,D and F levels separately is necessary because at 
the lower electron densities considered radiative de -excitation rates can 
become larger than collision rates between the n = 4 levels. At higher electron 
densities the collision rates between the n = 4 terms dominate all other rates 
in or out of these terms and the relative populations are Boltzmann. 


B. Population Equations 

The rate equation describing the population of the bound or continuum 
state i is 





) = 0 


II. 1 


5 

where Ry^i an ^ ®Yij are "total transition rates/cm by process y from state 
j to i and from i to j respectively. The sum S represents a sum over discrete 
states and an integration over continuum states. We shall assume the particle 
translational distribution"functions to be Maxwellian and the external 
continuum radiation field to be Planckian, in which case it is possible to 
integrate over the continuum and replace it by one additional term in the discrete 
sum. The atomic transition processes are radiative excitation and ionization, 
collisional excitation and ionization by atoms and electrons and their inverses. 
Because of their higher velocities, electron collisions generally dominate the 
collisonal rates and hence only electron inelastic rates will be considered. 
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TABLE II. 1 

He I, II Energy Levels 


He I 
3 

- i = 

1 

Energy (ev) 

Wave Nos. 

1 

Is 2 

x s 

0 

0 

2 

ls2s 

3 s 

19.821 

159850 

3 

ls2s 

1 s 

20.618 

166272 

4 

ls2p 

3 P° 

20.966 

I6908I 

5 

ls2p 

lpO 

21.220 

171129 

6 

ls3s 

3 s 

22. 721 

183231 

7 

ls3s 

1 s 

22.923 

184859 

8 

ls3p 

3 P° 

23.009 

185559 

9 

ls3d 

3 d 

23.076 

186096 

10 

ls3d 

■4) 

23.076 

186099 

11 

ls3p 

X P° 

23.089 

186204 

12 

ls4s 

3 s 

23.596 

190292 

13 

ls4s 

x s 

23.676 

190935 

14 

ls4p 

3 p 0 

23.710 

191211 

15 

ls4d 

3 d 

23.738 

191439 

16 

ls4d 

4d 

23.739 

191441 

17 

ls4f 

3 f° 

23.739 

191447 

18 

ls4f 

^F 0 

23.739 

191447 

19 

ls4p 

¥> 

23.744 

191487 

He II 

i = 

2 



1 

Is 

2 s 

0 

0 

2 

2s 

2 s 

40.8099 

329179. 57 

3 

2p 

2 P° 

40.8091 

329182.02 

4 

3s 

2 s 

48.3662 

390140.76 

5 

3p 

2 p 0 

48.3664 

390141.49 

6 

3d 

2 d 

48.3665 

390142.64 

7 

4s 

2 s 

51.0113 

411476.98 

8 

4p 

2 P° 

51.0114 

411477.28 

9 

4d 

% 

51.0115 

411477. 77 

10 

4f 

2 f° 

51.0117 

411477.95 


St.Wt. 

1 

3 

1 

9 

3 

3 

1 

9 

15 

5 

3 

3 

1 

9 

15 

5 

27 

7 

3 


2 

2 

6 

2 

6 

10 

2 

6 

10 


14 
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TABLE II. 2 
He I, II Lines 


He I 


Upper 

Level 

Lover 

Level 

Notation 

X( A) 

A 

(10^/ sec) 

f 


2 

1 

10830 

.1022 

-5391 

5 

1 

2 UV 

584.4 

17-99 

.2762 

5 

3 


20582 

. 01976 

.3764 

6 

4 

10 

7065 

.278 

.0693 

7 

5 

45 

7281 

.181 

.0480 

8 

2 

2 

3889 

.09478 

.06446 

8 

6 


4 . 30+4 

.0108 

.896 

9 

4 

11 

5876 

.706 

.609 

9 

8 


1.86+5 

1.28-4 

.111 

10 

5 

' 46 

6678 

.638 

.711 

11 

1 

3 UV 

537-1 

5.66 

.0734 

11 

3 

4 

5016 

_. i338 

.1514 

11 

7 


7-43+4 

. 00253 

.629 

12 

4 

12 

4713 

.106 

.0118 

12 

8 


21120 

.0652 

.145 

13 

5 

47 

5049 

.0655 

.00834 

13 

11 


21132 

.0459 

.103 

14- 

2 

3 

3188 

.0505 

.0231 

14 

6 


12538 

.00608 

.0429 

14 

9 


195^3 

.00597 

.0205 

14 

12 


I.09+5 

.0505 

1.21 

15 

4 

14 

4472 

.251 

.125 

15 

8 


17002 

.0668 

.482 

15 

14 


4.39+5 

4.16-5 

.200 

16 

5 

48 

4922 

.202 

.122 

16 

11 


19089 

.0711 

.647 

17 

9 


18688 

.139 

1.02 

17 

15 


1.43+7 

6.01-10 

.0033 

18 

10 


18699 

.138 

1.01 

18 

16 


1.67+7 

4.34-10 

.00253 

19 

l 

4 UV 

522.2 

2.46 

.030 

19 

3 

5 

3965 

.0717 

.0507 

19 

7 


15088 

.0137 

.140 
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TABLE II. 2 (Continued) 


He I 


Upper 

Level 

Lower 

Level 

Notation X( K ) 

A 

(10^/sec) 

f 

19 

10 

18560 

.00277 

.00858 

19 

13 

1.81+5 

5.78-4 

.853 

19 

16 

2.17+6 

5.65-7 

.024 


He II 


3 

1 

303,80 

100. 

.4162 

4 

3 

1640.5 

1.01 

.01359 

5 

l 

256.3 

26.8 

.07910 

5 

2 

1640.4 

3.59 

.4349 

6 

3 

1640.4 

10.3 

.6958 

7 

3 

1215,1 

.413 

3.045-3 

7 

5 

4687.0 

.294 

.03225 

8 

1 

243.03 

10.9 

.02899 

8 

2 

1215. 1 

'1.55 

.1028 

8 

4 

4686.8 

.491 

.4847 

8 

6 

4687.2 

.056 

.01099 

9 

3 

1215.1 

3.3o 

.1218 

9 

5 

4686.9 

1.13 

.6183 

10 

6 

4687.1 

2.21 

1.018 
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We may write the total rate from i to j more specifically as. 


R . . = n.P. . = n.(A. . + C. .), 


II.2 


where A., and C. . are the radiative and collisonal transition rates/particle 
ij 

from the ith to the jth state. In the statistically steady state, equation 
II. 1 reduces to 


S 


(D /ji 


- n iV 


= n /ji ' 0J 




II. 5 


The system of equations represented by II . 3 can easily be solved for the 

populations n. when all the P. .'s are known. This is the case when the gas 
i ij 

is totally optically thin (in all lines), since the external radiation field 
is specified and the internal radiation field does not produce a significant 
upward transition rate. We can characterize the system of linear equations 
(II. 3 ) by a matrix whose coefficients a^ are equal to P . The diagonal 
elements are the P. . . 

li 


The general solution of equations (II. 3 ) is given by White (1961) 


n. 

1 




II. 4 


i 

where P is the co-factor of the coefficient of n . in the mth equation 
(i.e. the matrix element P^) and N is the total number of He particles/cm . 

We denote the discrete level corresponding to the continuum, i.e. an ionization 
or recombination, by c. We may then write the rate equation for level i of 
ion q as. 


p i/ + - Yn n j= 0 

J J J J J 

where j and j 1 refer to other bound levels in the qth ion and in other stages 

of ionization, respectively, P^. is the total excitation or de-excitation rate 

q c 

in ion q from i-j; P3 . , and P.T, . are the total ionization rates from i in 

J 1 fl 

q to j ’ in q + 1 and to i in q from j 1 in q - 1, respectively; P . and 

C1J 
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p?^ + ^ are the total 
in q from j' in q+1, 


recombination rates from i in q to j' in q -1 and to i 
respectively. Specifically we may write the P's as follows: 


bound-bound 

i > 0 
i < 0 


= A?. + Y^.B^.B +n0?. 


ij 


iJ 


ij ij r e ij 


pj, = Y^.B^.B + n C& . 
ij ij ij r e ij 

bound -free (from initial state q, i) 


II. 6 


ionization 

pi 

icj' 


i j' 



II. T 

recombination 

pi 

cij* 

■ ViJ- + 

n e W PiJ' 4 

n^f 

e 

q 

ij r 


bound -free 

(to initial state q 

1) 




ionization 

q-l,c 

■ w5 Ji 4 




II. 8 

recombination 

pC , q+1 
j'i 

* n e a r i 4 

' n e WP J'i 4 

2 

n 

e 

cf i +1 
j’i 



where is the electron density/cm ; is the Planck function at temperature 

T ; W is the dilution factor; Yj*\ is a free parameter. A? 1 ,, B?-. are the Einstein 

r - q ij lJ lJ a - a 

transition probabilities; A ... is the photoionization rate/ion; U.., fl . . are 


the collisions! transition_rates/electron per ion for bound-bound and bound-free 

processes , respectively; is the collisional transition rate/electron per 

ion for free -bound recombination; and C$,, f}1. are the recombination, and 

ij *ij ’ 

stimulated recombination coefficients, respectively. 


The basic quantities required for the evaluation of the radiative rates 
are the oscillator strengths and photoionization cross sections; for the colli- 
sional rates the excitation and ionization cross sections. The inverse cross 
sections and Einstein coefficients can be obtained from the usual detailed balance 
relations. The quantities U, U, U and a generally have the form 


,f ^(v)/(v) dv, 


XI. 9 



12 


where v is the electron velocity, Q,(v) the cross section and /(v) the electron 
translational distribution function. A may be written in terms of the photo- 
ionization cross section a( v) as 

4ir J B(v) dv. II. 10 

^ is expressed similarly in terms of the stimulated recombination cross section 
b(v). 


Solutions can be obtained for any arbitrary line radiation field by 
suitable choice of the parameter Y which is closely related to the net radia- 
tive bracket (NRB). Y varies between 0 and 1 being 0 for a thin layer (no 
external radiation) and 1 for a very thick layer where the line radiation field 
is Planckian. For a thin line with external photospheric radiation Y = l/2 B(T e ) . 
Various continuum radiation fields can be considered by suitable choices of 
A. 



13 


C. Reaction Rates 

Oscillator strengths for all of the allowed transitions included in our 

model (Section II. A) were found in the KBS compilation by Wiese et al. (19 66). 

Photo ionization rates were obtained from Hartree-Fock calculations of Stewart 

and Webb (1963) for the ground state of He I and from calculations by Peach (1967) 
13 1 3 

for the n' S, n' S, n P and n P levels. For all other levels cross sections 
were calculated using the quantum-defect method of Burgess and Seaton (i960). 

Collisional ionization rates were obtained from measurements of 
Englander -Golden and Rapp quoted by Kieffer and Dunn (1966) for 1 S, from measure- 
ments of Long (1967) for 2^ S and from calculations of Bolder, et al. (196$ 
for the He II ground state. For other levels the ionization rates were taken from 
Allen (1961). 

The collisional excitation cross sections are crucial to the solution 
of the steady state equations. For this reason we have made an extensive 
tabulation and study of rates from many sources. A paper entitled "Electron 
Impact Excitation Rates for Helium" describing these rates was published in 

the December 1972 issue of Solar Physics (Benson and Kulander, 1972). Excitation 
rates were calculated from modt available cross section data, and fitted to the 
empirical formula 

n = AT n exp (- a X Q ) 

where X = E /kT ; A, n and Q£ are constants. For He I the temperature range 

00 CL 4 

considered was 4000-50, OOOTC and for He II, 10 - 10 K. Rates between all 

levels of the model of Section II. A were calculated. The inverse rates for both 
radiative and collisional transitions were calculated from standard equilibrium 
relationships . 

We shall discuss briefly some of our conclusions concerning the electron 
impact rates beginning with the He I rates. 



lb 


Generally both forbidden and allowed rates are in better agreement 
for higher temperatures than for lower temperatures. The allowed rates are 
in much better agreement than the forbidden rates. The allowed rates from 
ground state or level 2 generally show differences as large as a factor of 
10 for low temperatures and as large as 5 for high temperatures. Between 
other levels (n ^ 3) the differences are only as much as factors of 3 at 
low temperatures and 50-100$ at high temperatures. 

In He I the forbidden rates generally differ by as much as a factor 

of 100 at low temperatures, a factor of 20-50 at higher temperatures. Differ- 

5 

ences as high as a factor of 10 are noted in a few cases, differences of 

3 

factors of 10 are not uncommon. If one does not consider the highest and 
the lowest rate value for each transition the differences are generally 
reduced to factors of about 10 at lower temperatures and 5 at higher tempera- 
tures. Forbidden transitions with lower level n = 2 show differences 2-3 
times less than those with n = 1. 

For transitions between the higher levels only the very approximate 
cross sections of Green (I966) and Allen (1963) are available for forbidden 
transitions and for allowed transitions those of Seaton (1962), Saraph (1 964) 
and Mihalas and Stone (1968). Comparing the two approximate results for 
forbidden transitions with the other measured and calculated values we find 
the approximate results generally lower by factors of 2-10. For higher levels 
the Green cross sections are generally higher than those of Allen. With re- 
gard to the Seaton and Saraph approximations, it is not conclusive which is 
more often closer to the other values. The Seaton approximation gives values 
generally higher than the Saraph approximation for lower temperatures while 
the reverse situation holds for higher temperatures. 

For He II there are relatively few cross sections available. The 
largest differences, being about a factor of 5> are much less than for He I. 
For both He I and He II experimentally determined cross sections give lower 
rates than the calculated values. 
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The actual collisional excitation rates chosen represent mean values 
of those given in the paper, taken from various sources listed there. To 
determine the sensitivity of the solution to these rates the statistical equi- 
librium equations were solved after increasing the individual values of 
the collision rates per electron Q between each state by a factor of 100. This 
was done to determine the sensitivity of the solution to the collision exci- 
tation rates. We chose the optically thin case with T = 20,000°K, n g = 10 10 
for this test. We illustrate in Figure II. 3 the effect on the population of 
the levels 1, 2, 4, 9, 11 and 20 of each perturbation in collision rate. Of 
course in many instances increasing a particular Q. . by 100 had little effect 
on the populations of the levels mentioned. In Figure II. 3 we show the 10 
transitions i-j which have the largest effect on each of the above levels. 

The ratio of the population after increasing by 100 to that before the 
perturbation is shown. We note that changes of 10 to 50 in these important 
populations result from an uncertainty of a factor of 100 in the rates. 

All the triplet levels of He I had about the same response to the change 
in Q. When the rate from the He I ground state (G.S.) to triplet levels increased 
the triplet populations increased. When the triplet to singlet rates increased 
the triplet populations decreased as would be expected. The He I G.S. population 
is very sensitive to the rate from the G.S. to the metastable 2 S level &nd 
2 level). The 1 level population was decreased by about 20 when 
fl(l S — 2 -») was increased by 100. At the same time the He II G.S. pupulation 

was increased. An increased rate from the triplets to singlets results in 

1 2 
increased 1 S 'population and a corresponding decrease in 1 S population. 

The populations of the upper levels of the 58k and 557 A lines are very sensitive 

to the rate from the 1 to 2 1 S level increasing as this rate increases. 
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D. Sample Solutions 

A code which we shall designate as code 1 has been developed to solve 
the steady state population equations given in Section II. B. This code is 
described in Section V. We discuss in this section some sample solutions of 
the steady state population equations. 

The gas is assumed irradiated over 2 k ster by a Blackbody spectrum 
at 6000°K representing the photospherlc radiation field. The gas is assumed 
to be optically thin for all lines and continuua unless otherwise specified. 

1. Ionization Equilibrium 

In this section we illustrate results for the ionization equilibrium. 

B.g.n-4 shows the ratio n^/rtg- with T & = 30,000°K at various values of n g . 

For an atomic model with only levels 1 and 20 the ratio n 1 / n 2o is approximately 

constant with n since both the collisional ionization and recombination rates 
e 

are proportional to n . The optically thin solution for two levels is shown. 
The results are generally within a factor of 2. For very low n g , all of the 
levels of He I except 5, 11 and 19 represent additional paths from level 1 
to level 20. This is so because the photoionization rates from these levels 
exceed almost all collisional rates. n-/n_,_ is hence lower at n =10 than 
the two level solution. For somewhat higher electron densities the collisional 
rates between the singlets and triplets exceed the photoionization rates. 

Hence, the recombination to intermediate levels results in conversion to the 
5, 11 or 19 level and thence to the 1 level by spontaneous emission. The 
rate from 20 to 1 is now enhanced and n^/ exceeds the two level solution. 

At still higher electron density the collisional ionization rates begin to be- 
come greater than either the spontaneous emission or collisional deexcitation 
rates. Now all intermediate levels represent paths from level 1 to 20 and 

n -i/hpn becomes lower than the two level solution. This is seen to be the case 

1 14 

at n = 10 . 

e 

We turn to the solution when the resonance lines are optically thick. 
The result for detailed balance (db) in the first resonance line is shown 
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by the curve labeled = 0 (p is the net radiative brackett). Also shown are 
results for db in the second and third resonance lines. Increasing the radiative 
excitation rate increases the effective ionization rate and the amount of n^Q by 
several orders of magnitude. We thus note the important result that the ioniza- 
tion equilibrium is strongly influenced by the optical thickness of the resonance 
lines. The LIE result is also shown and is labeled LTE. 

2 . Level Populations 

Many sample solutions have been obtained for the level populations or non- 

equilibrium parameters, b, where b = n/n 

' equil 

We have completed a parametric study of the effect of optical thickness in the 

resonance lines and continuua upon the level populations. Level populations were 
obtained for various physically meaningful combinations of resonance lines and 
continuua being optically thick or thin. That is to say each resonance line or 
continuua was assumed to have a net radiative brackett of either 0 (completely 
thick) or 1 (thin). The results cover the temperature range T = 10,000 to 
50,000°K and electron density range n g = 10 10 -10 l4 cm" 5 . Due to space limita- 
tions we cannot present results for all 30 levels of the model. Hence levels 

I, 5, 9, 20 and 27 were chosen to illustrate the results. Referring to section 

II. A we see that levels 1 and 20 are the ground states of He I gr>H ji ? 

level 5 the upper level of the 59^-A line, level 9 the upper level of the D3 line 

and level 27 the upper level of the 1+686A . line of He II. Figures II. 5-29 show 

results for each of these 5 levels for electron temperatures lo\ 2 x 10^ 

3 x 10 , b x 10 and 5 x 10 K. In each figure the ratio of the actual popu- 

lation to the equilibrium population, b, is given as a function of n . 

e 

At each temperature calculations were made for combinations of net radiative 
brackets corresponding to layers of varying total thickness. There are 5 
resonance lines in the model for both He I and He II. These 6 lines together -Kith 
the Lyman continuua for each ion are allowed to become optically thick in 
our calculation. Thus, there are eight lines and continuua which can be 
optically thick or thin depending upon the physical thickness of the layer. 

Each figure shows a completely optically thin solution labeled T and 8 
other solutions labeled 1-8 which correspond to the combinations of optical 
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Figure II-6 


Te = 10,000° K 
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Figure II -10 


> Te - 20,00(fK 

A - - . - - 7 _ 
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Figure II-lU b . 
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thicknesses given in Table II. 3. On this table (l) refers to the 5841 line, (2) 
refers to the 537A line, (3) refers to the 522A line and (4) refers to 
the Lyman continuum of He I. Similarly the numbers 5-8 refer to the 304, 

256, 2 4^ A lines and Lyman continuum of He II. When a number appears in 

Table II. 3 the corresponding transition has been assumed optically thick 
(i.e. in radiative detailed balance) in obtaining rate coefficients for the 
level population solutions. Transitions not appearing in the table are 
assumed optically thin. It is noted that progressing from case 1 to 8 corres- 
ponds in general to the layer becoming thicker. Case 8 always corresponds to 
all 8 lines and continuua becoming optically thick. For example, with 
T e - 50,000°K case 4 represents the NRB = 0 in the first 5 resonance lines of 
He I and the first resonance line of He II. 

We shall discuss briefly the solutions presented in Figures II.5-3Q. For 
4o 

T g = 10 Kmost of the He is He I. The He I ground state population is not 

affected by optical thickness except slightly at n = 10 10 . Similarily Figure II. 6 

0 e 

shows that once the 584 A line becomes thick the radiation field in other lines 

and continuua do not affect it. The D3 line upper level population shown in Fi- 

gure II . 7 is a strong function of the 584 A line and He I Lyman continuum optical 

o 

depths and to a lesser extent dependent on the 537 and 522 A lines. Increasing 
the singlet populations hence increases the D3 emission. The optical thicknesses 
of the He II lines and continuua do not affect the D3 line. From Figure II. 8 it is 
seal that the He II ground state exhibits a similar behavior. At 20,000°K the 
He I and He II densities are very roughly the same. The thin solution shows 
that He I is greatly overpopulated and He II slightly underpopulated with respect 
to LTE. Increasing the optical depth in the He I resonance lines and continuua 
now increases the ionization by absorption of photospheric radiation from upper 
resonance line levels and thus increases the He II density and decreases the 
He I density. The D3 line upper level population increases with increasing 
population of the singlet levels but decreases when the He I Lyman continuum 
becomes thick due to decreasing total He I density. As before the He II resonance 
line and continuum radiation fields do not affect the He I level populations. 

At 30, 000° K both He I and II ground states are overpopulated at the expense 
of He III. The He II ground state population is sensitive only to the Lyman 
continuum of He II while the He I population depends upon both continuua and the 
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Table II -3 

OPTICALLY THICK LIMES CHOSEN 


Case Number 

1 2 5 4 b 6 7 8 



He I resonance lines. The D3 line is sensitive only to the optical thickness in 
the two Lyman continuua. The level 27 population depends upon the He II resonance 
lines and continuua only. At T g = 40, OOCP K all the levels have a strong dependence 
on the optical thickness of the He II Lyman continuum. When this continuum is thick 
the population of He III is greatly increased while those of He I and II are decreased, 
then decreasing all of the level populations. For levels 9 and 20 this continuum 
has the largest influence. Levels 1 and 5 are also dependent on the He I resonance 
line and continuum optical depths. At 50,000° K the level populations have the 
same behavior as at JiOjOOCf’K. 




k7 


3. Optical Thickness and Line Intensity- 

Figures 11.5-29 have illustrated the effect of optical thickness in 
the resonance lines and continuua on certain level populations. We now seek 
to determine the approximate optical depths of a flare layer. For purposes 
of illustration we consider a 1000 km thick layer. For each of the same cases 
given in Table II. 3 we have obtained the line center optical thickness in the 
58k, 304, D3, 10830A lines and at the threshold of the Lyman continuua of 
He I and II. Optical depths are not shown for the 537, 522, 256, 243 A lines 
because they are always simply a constant fraction of the 584 and 304A 
line optical thicknesses, namely t(537) = .27 t(584), t(522) = ,11t(584), 

t(256) = .19 t(304), and t(243) = .070 t(304). 

The optical thicknesses are given in Tables II. 4 - II. 8 The He 
density is assumed to be one tenth of the electron density. The number 
following each entry is the power of 10 by which the entry is multiplied. 

We note that the 10830 and D3 lines can become thick for high electron den- 
sities even at 50,000°K. These lines do not become thick at 10,000°K. There 
are many cases in which a number of lines and continuua are optically thick. 
This does not mean however that simultaneous transport equations must be solved 
for these lines and continuua. Which line and continuum radiation fields must 
be obtained simultaneously depends upon the level population being sought as 
well as the temperature, density and layer optical thickness. To determine the 
effect of various lines and continuua on levels 1, 5, 9, 20 and 27 we can refer 
to Figures II. 5 - 29 . For example. Table II. 9 shows the approximate maximum 
percentage error encountered in the I>3 line upper level population by solving 
only for the lines and continuua given in each box. The effect on a^y level 
can be obtained from the output from Code 1. Figures II . 5 - 29 illustrate the 
effect of varying the NEB between 0 and 1. In practice when the layer becomes 
optically thick it does not usually become so thick that the HRB = 0. In order 
for the line radiation field to become saturated we require 

> (e + Tl)" 1 , 
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Table II -4 

OPTICAL THICKNESS - 1000 KM LAYER 

T = 10,000 °K 
e 


Case 


TRANSITION 

n 

- £ 

T 

1 

2 

3 

4 

5 

6 

7 

8 


10 10 










584 


3. 7+^ 

- 

5S 

a 

1. 

= 

a 

= 

= 

504 


7.6 


= 

= 

4.1 

= 

* 

= 

- 

5876 


2.8-8 

1.1-5 

2.1-5 

2.7-5 

8.0-4 

= 

= 

= 

= 

10850 


2.6-7 

1.0-4 

1.9-4 

2.5-4 

7.4-3 

= 

= 

= 

= 

504 


1.5-2 

9-6 



8.6+2 

= 

= 

ss 

= 

228 

10 11 

1.4-6 

9-3-4 

1.8-3 

2.4-3 

8.3-2 

x: 

rs 


ss 

584 


3.7+5 

= 

- 

= 

= 

= 

= 

= 

= 

504 


76 

= 

= 

= 

= 

= 

as 


- 

5876 


1.3-6 

1.7-4 

2.2-4 

3.7-4 

6.7-3 

= 

3* 

= 

= 

10830 


1.2-5 

1.2-3 

2.1-3 

3-4-3 

6.2-2 

= 

= 

= 

sc 

304 


.046 




1.6+3 

- 

S3 

S3 

= 

228 

10 12 

4.4-6 

1.0-3 

1.9-3 

3-1-3 

.15 



= 

= 

584 


3-7+6 

= 

= 

= 

= 

= 

= 

= 

- 

504 


7.6+2 

= 

= 


= 

a 


SS 

a 

5876 


1.5-5 

1.1-3 

2.4-3 

4.7-3 

1.5-2 

= 

- 

= 

* 

10830 


1.2-4 

8.6-3 

2.0-2 

3.7-2 

0.12 

rs 


- 

= 

504 


0.10 




1.3+3 

a 

= 

ss 


228 

10 13 

9.3-6 

1.0-3 

2.8-3 

6.0-3 

1.2-1 

= • 

= 

sc 

— 

584 


3.7+7 


= 

= 

= 

= 

= 

= 


504 


7.6+3 

= 

= 

= 

S5 


= 

= 

= 

5876 


2.1-4 

1.0-2 

3.1-2 

4.1-2 

6.2-2 

. = 

= . 

= 

= 

10330 


1.4-3 

4.7-2 

0.13 

0.17 

0.26 

=3 

= 

= 

= 

304 


0.56 


1.3+22 

2.5+2 

1.1+3 

- 

= 

= 


228 

lO 1 ' 1 

5.4-5 

2.7-3 

1.3-2 

2.4-2 

1.0-1 



= 

= 

584 


5,718 

~ 

= 

- 

= 

- 

= 

- 

= 

504 


7-6+4 

= 

= 

= 

= 

a 

= 

= 

- 

5876 


5-6-3 

0.15 

0.24 

0.26 

0.29 

= 

= 

= 

=s 

10830 


1.4-2 

0.37 

0.54 

0.59 

0.66 


- 

= 

= 

304 


4.7 

1.7+2 

5.0+2 

6.1+2 

9.3+2 

- 


= 

= 

228 


4.5-4 

1.6-2 

4.8-2 

5.9-2 

9.0-2 

a 


= 

= 
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Table II -5 

OPTICAL THICKNESS - 1000 KM LAYER 

k rs 

T = 2 x 10 K 


Case 


TRANSITION 

n 

e 

T 

1 

2 

3 

4 

5 

6 

7 

8 


10 10 










584 


1.7+3 

3.4 

0.92 

0.52 

= 

= 

= 

7.0-3 

= 

504 


0.56 

7.2-4 

1.9-4 

1.1-4 

= 

- 


1.5-6 

= 

5876 


1.7-3 

1.8-3 

1.7-3 

= 

= 

= 

= 

9.6-4 


10850 


1.6-2 

1.7-2 

1.6-2 


= 

= 

= 

9.3-3 


504 


1.5+3 

2.8+3 

= 

- 

= 

= 

- 

2.0+3 

= 

228 

10 11 

0.10 

0.19 

St 


= 

~ 

= 

0.13 

= 

584 


2.0+4 

2.2+2 

61. 

26. 

= 


ss 

1.0 

= 

504 


5-7 

6.3-2 

1.7-2 

7.4-3 


= 


2.9-4 


5876 


4.5-2 

0.11 

= 

= 

= 


= 

3.3-2 

= 

10850 


0.58 

1.0 

0.86 

0.93 

= 

' = 


2.9-4 

= 

504 


5.1+3 

2.1+4 

- 

= 


= 

“ 

= 

- 

228 

10 12 

0.45 

1.8 

— 


= 

— 

~ 

” 

— 

584 


2.3+5 

1.0+4 

2.0+3 

5.6+2 


= 


100. 

= 

504 


65. 

2.8 

0.56 

0.16 



= 

2.8-2 


5876 


0.90 

4.3 

4.0 

3.1 

= 

s 

= 

1.0 

= 

10850 


5.3 

24. 

22. 

17. 


= 

= 

= 


504 


3.3+4 

2.0+5 

= 


= 

= 

= 

= 


228 

10 15 

2.7 

16. 

— 

' 

= 

— 


= 


584 


2 . 0+6 

1.2+5 

1.7+4 

8.5+3 


= 

= 

5.2+3 

= 

504 


6.1+2 

37. 

5.2 

2.6 

= 

. = 

= 

1.6 


5876 


18. 

63. 

48. 

38. 

= 

' = 

= 

28. 

= 

10850 


37.. 

1.5+2 

90. 

69. 


= 

r= 

52. 

= 

504 


4.0+5 

1.9+6 


- 7 s 

= 

= 

= 

= 

= 

228 

io 1 ^ 

36. 

H 

0 




— 



— *■ 

584 


2.0+7 

1.2+6 

3.8+5 

3.0+5 

s 

- 


2.6+5 

= 

504 


5.8+3 

350. 

no. 

87 . 

:= 


- 

75. 


'iflyo 



7.7*2 

5.4*:’ 

5 . 1 +:' . 


=: 


4.8+2 

= 




v'.i').' 

5.8+2 

5 . 4+2 

-- 

- 


5.1+2 

- 

soh 



l.'M-y 

= 

= 

= 

= 


= 


228 



1.7+5 

= 

= 

= 
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Table II -6 

OPTICAL THICKNESS - 1000 KM LAYER 

T = 3 x 10 °K 
e 


Case 


TRANSITION 

n 

e 

T 

1 

2 

3 

4 

5 

6 

7 

8 


10 10 










584 


25 . 

= 

22 . 

21. 

'2.5-2 

5 . 2-3 

2 . 7-3 

1 . 7-5 

1.6-5 

504 


8 . 4-3 



8.2-3 

1.0-6 

2.0-7 

1.0-7 

5.9-6 

6.0-11 

5876 


1 . 3-5 

= 

= 

1.2-3 

7.7-^ 

7 . 0-4 

= 

4.3-6 

3.2-6 

10830 


1.2-2 


- 

1.1-2 

7.2-3 

6 . 5-3 

6 . 4-3 

4.0-5 

3 . 0-5 

304 


1 . 6+3 

- 

= 

= 


= 


10. 

= 

228 

10 11 

0.18 


- 

0.17 

= 

= 

= 

1 . 1-3 

— 

‘ 584 


6.3+2 

= 

= 

- 

2.0 

0.43 

0.15 

8.8-3 

2 . 6-4 

504 


0.19 

:= 

= 

= 

7 . 2-4 

1 . 5-4 

5 . 0-5 

3.2-6 

1.0-7 

5876 


6.8-2 

= 


= 

5.2-2 

4.5-2 

= 

2 . 9-3 

8.2-5 

10830 


0.57 

= 

= 

= 

0.43 

0.37 

= 

2 . 4-2 

7 . 0-5 

304 


1 . 6+4 


- 

= 

- 

= 

= 

1.0+3 

= 

228 

10 12 

1.7 

— 

— 

~ 

= 

— 

— 

0.10 

= 

584 


7 - 5+5 



= 

62. 

9.6 

2.5 

1.0 

0.13 

504 


2.6 

= 


= 

.22 

3 . 4-3 

8 . 9-4 

3 . 6-4 

4 . 6-5 

5876 


1.5 



- 

1.4 

= 

1.1 

0.42 

0.15 

10830 


6.8 

= 

- 

= 

6.2 

5.9 

4.7 

1.8 

0.63 

304 


1.6+5 


= 

= 

= 

: = 

= 

6 . 4+4 

= 

228 

10 15 

18. 



= 

— 

= 

— 

7.0 


584 


4 . 8+4 

= 

= 

- 

6.5+2 

76. 

38 . 

32 . 

17. 

504 


17 . 


=Z 

= 

0.23 

2.7-2 

1 . 4-2 

1.2-2 

6. 4-3 

5876 


22. 

= 

a; 


18. 

15 . 

12. 

10. 

7.3 

10830 


52 . 


= 

= 

'25. 

19 . 

15 . 

13 - 

10. 

304 


1.6+6 

= 

- 

' = 

= 

= 

- 

1 . 4+6 


228 

10 lU 

1.8+2 

2 = 

- 

= 

= 

= 

= 

1.5+2 

= 

584 


5 . 1+5 

50+5 

3 . 7+5 

= 

5.^3 

1 . 8+3 

1 . 4+3 

= 

1.1+3 

504 


1.2+2 


10+2 

= 

2.3 

O.69 

0.48 

= 

0.38 

5876 


2.0+2 


2 . 4+2 

= 

1.9+2 

1 . 4+2 

1 . 3+2 

= 

1.2+2 

10830 


2.6+2 

= 

2 . 1+2 

= 

1.7+2 

1 . 2+2 

1 . 1+2 

= 

1 . 0*2 

304 


1 . 5+7 


- 

= 

1.6+7 

= 

s 


= 

228 


1 . 8+3 

= 

- 

= 

= 



— 

- 
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Table II -7 

OPTICAL THICKNESS - 1000 KM LAYER 

T = 4 x 10^ 0 K 
e 


Case 


TRANSITION 

n e 

T 

1 

2 

3 

4 

5 

6 

7 

8 


O 

O 

i — 1 










584 


1.9 


1.5 

0.37 

1.9-4 

1.0-7 

1.7-8 

8.5-9 

5.5-11 

504 


7.6-4 

- 

4.5-4 

1.5-4 

7.6-8 

4.0-11 

6.8-12 

3.4-12 

2.2-14 

5876 


8.3-4 

- 

6.4-4 

1.6-4 

1.6-8 

1.4-8 

1.2-8 

= 

6.6-9 

10830 


7-7-3 


5.9-3 

1.5-3 

7-7-7 

1.3-7 

1.0-7 


6.0-8 

304 


1.5+3 

= 

1.2+3 

2.9+2 

2.9-2 

as 

= 

= 

= 

228 

10 11 

O.18 


0.11 

.035 

3.5-6 

— 


' 

— 

584 


4o. 

38. 

33. 

15. 

6.7-3 

4.3-5 

8.0-6 

2.7-7 

4.7-9 

504 


1.7-2 

ss 

1.5-2 

7.0-3 

2.8-6 

1.8-8 

3 . 4-9 

l.l-io 

2.0-12 

5876 


O.32 

0.31 

0.27 

0.13 

5.6-6 


3.8-6 

= 

1.2-6 

10830 


0.26 

0.25 

0.21 

0.11 

4.5-5 

= 

3 . 9-5 

= 

1.3-5 

304 


1.4+4 

- 

1.3+4 

7.6+3 

2.8 

= 

- 

= 

= 

228 

10 12 

1.8 

— 

1.6 

0.81 

3.6-4 

— 

— 


— 

584 


4.2+2 

4.1+2 

3.7+2 

2.5+2 

0.82 

8.1-3 

1.1-3 

2.5-4 

3.3-5 

504 


0.17 

- 

0.15 

0.11 

3.4-4 

3.4-6 

4.7-7 

1.2-7 

1.4-8 

5876 


0.69 

0.67 

0.61 

o.4i 

1.4-3 

= 

* 

1.1-3 

4.0-4 

10830 


2.6 

2.5 

2.3 

1.5 

5.0-3 

= 

4.9-3 

4.0-3 

1.4-3 

304 


1. 5+5 

1.4+5 

1.3+5 

8.7+4 

2.9+2 

= 

= 

= 


228 

10» 

18. 

17. 

16. 

11. 

3.6-2 


=s 

= 

= 

584 


2.7+3 

2.5+3 

2.3+3 

1.6+3 


0.80 

0.10 

4.5-2 

2.0-2 

504 


1.0 

0.92 

0.84 

0.68 

2.2-2 

2.1-4 

3.6-5 

1.8-5 

7*7-6 

5876 


9.1 

8.7 

8.0 

5.6 

0.19 

0.16 

0.13 

0.10 

8.0-2 

10830 


11. 

10. 

9.2 

6.5 

0.22 

0.19 

0.l4 

0.12 

8.5-2 

304 


1.5+6 

1.4+6 

1.2+6 

8.8+5 

3.0+4 

= 


= 

= 

228 

10 lU 

180. 

170. 

160. 

100. 

3.6 

= 

= 

= 

= 

584 


2.1+4 

2.0+4 

1.7+4 

1.2+4 

3.1+3 




9.0 

504 


8.4 

7.8 

6.4 

3.0 

0.77 

1-3-2 

4.0-3 

3-2-3 

2.2-3 

5876 


90 . 

05. 

83. 

58. 

15 . 

12. 

10. 

9.0 

8.3 

1 0830 


76 . 

72 . 

62. 

45. 

11. 

10. 

7.0 

6.6 

6.0 

30J1 


1 . 4+7 

- 

t .2+7 

9.0+6 

2.3+6 

s 


= 

= 

228 


1 .8+3 

1.7+5 

1 . 4+3 

1.2+3 

3-0+2 

= 

=s 

= 

— 
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Table II -8 

OPTICAL THICKNESS - 1000 KM LAYER 

L. 

T = 5 x 10 °K 
e 


Case 


TRANSITION 

n e 

T 

1 

2 

3 

4 

5 

6 

7 

8 


10 10 










584 


0.35 

0.21 

4.2-2 

3.5-3 

2.1-7 

- 

- 

- 

- 

504 


1.5-4 

9.0-5 

1.8-5 

1.5-6 

9.0-11 

- 

- 

- 

- 

5876 


5.4-4 

3.4-4 

6.5-5 

5.3-6 

3.2-10 

- 

- 

- 

- 

10850 


4.9-5 

3.0-3 

6.0-4 

5.0-5 

3.0-9 


- 

- 

- 

304 


1.3+5 

7.8+2 

1.6+2 

13. 

7.8-4 

- 

= 

= 

= 

228 

1 — ! 

1 — 1 

O 

H 

0.18 

0.11 

2.2-2 

1.8-3 

1.1-7 

- 


= 

= 

584 


6.6 

4.0 

1.3 

0.26 

4.0-5 

- 

- 

- 

- 

504 


3.0-3 

1.8-3 

5.9-4 

1.2-4 

1.8-8 

- 

- 

- 

• 

5876 


1.8-2 

1.1-2 

3.6-3 

7.2-4 

1.1-7 


- 

- 

- 

10830 


0.l4 

8.4-2 

2.8-2 

5.6-3 

8.4-7 

- 

- 

- 

- 

304 


1.3+4 

7.8+3 

2.6+3 

5.2+2 

7.8-2 

= 


= 

- 

228 

10 12 

1.8 

l.l 

0.36 

7.2-2 

1.1-5 


- 

sr 

= 

584 


62. 

37- 

12. 

3.1 

3.7-3 

- 

- 

- 

- 

504 


3.0-2 

1.8-2 

5.8-3 

1.5-3 

1.8-6 

- 

- 

- 

- 

5876 


0.37 

0.22 

7.4-2 

1.8-2 

2.2-5 

- 

- 

- 

- 

10830 


1.2 

0.75 

0.24 

6.0-2 

7.2-5 

- 

- 

- 

- 

304 


1.3+5 

7.8+4 

2.6+4 

6.5+3 

7-8 

- 

= 

s 

= 

228 

10 15 

18. 

11. 

3.6 

0.90 

1.1-3 

ZZ 

zz 

— 


584 


4.1+2 

2.4+2 

82. 

20. 

O.25 

4.0-3 

4.5-4 

2.3-4 

9.0-5 

504 


0.18 

0.11 

3.6-2 

8.8-3 

1.1-4 

1.8-6 

2.0-7 

1.0-7 

4.0-8 

5876 


4.8 

2.9 

0.96 

0.24 

2.9-3 

2.6-3 

2.2-3 

1.7-3 

1.3-3 

10830 


4.9 

2.9 

1.0 

0.25 

2.9-3 

2.6-3 

2.1-3 

1.7-3 

1.3-3 

304 


1.3+6 

7-&+5 

2.6+5 

6. 5+4 

7.&+2 

= 

= 

= 

= 

228 

10 14 

180. 

no. 

36. 

8.8 

0.11 


- 

= 

— 

584 


3.0+3 

1.9+3 

5.0+2 160. 

19. 

0.35 

0.11 

8.5-2 

5.9-2 

504 


1.4 

0.84 

0.22 

7.0-2 

8.9-3 

1.6-4 

5.1-5 

4.0-5 

2.8-5 

5876 


47. 

30. 

7.8 

2.4 

0.29 

O.25 

0.20 

= 


10830 


33. 

21. 

5.6 

1.7 

0.21 

0.18 

0.14 

= 

- 

304 


1. 5+7 

7.8+6 

2.1+6 

6.5+5 

7.8+4 

= 

= 

= 

= 



1.8-0 

•1.1+5 

i?.y+L? 

9.0+1 

11. 

= 

- 


XX 




Table II -9 ERROR IN D 3 Line Intensity 


fic X “ D 3 Line 


io* V 


iXiO 4 °K 


3 \\ o H ° K 


*» o 

4 X « 


10 * ° K 



* 6 t A. I 13 *4 T 6 

fi. Lime | ST 1 N ° 3*° L M~c I s ~^ 1 **■ 
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TABLE II. 10 

CHARACTERISTIC VALUES OF e, T\, C 


T (°K) 
e 

u e , 

G 


i 

Line 

(A) 

Case 

# 

10 4 

10 10 

4.1-9 

2.3-4 

2.9-18 

584 

T 



tt 

Tl 

5.1-14 

It 

3 



1.2-8 

1.2-8 

1.8-28 

304 

T 



ft 

II 

tt 

It 

3 


10 14 

4.0-5 

8.0-3 

4.0-14 

584 

T 



II 

II 

8.3-13 

ft 

3 



1.2-4 

4.8-5 

1.8-24 

304 

T 



ft 

If 

If 

II 

3 

4. 

2 x 10 

10 10 

1.1-8 

2.3-4 

6.1-13 

584 

T 



II 

If 

3.8-3 

tt 

3 



8.2-9 

1.2-8 

2.5-18 

304 

T 



It 

fl 

tl 

tt 

3 


10 14 

1.1-4 

2.2-2 

6.4-9 

584 

T 



ft 

tl 

4.9-7 

tt 

3 



8.2-5 

5.8-5 

2.5-14 

304 

T 



tl 

T? 

II 

n 

3 

3 x 10 4 

10 10 

1.9-8 

2.3-4 

4.3-11 

584 

T 



Tl 

91 

II 

It 

3 



6.7-9 

1.4-8 

6.1-15 

304 

T 


10 14 

II 

tl 

4.2-9 

tt 

3 


6.7-5 

1.2-4 

5.8-11 

584 

T 



tl 

M 

1.4-8 

tt 

3 



I.9-4 

3.4-2 

4.0-7 

304 

T 



n 

II 

It 

tt 

3 

4 . 

4 x 10 

10 10 

2.8-8 

2.3-4 

3.9-10 

584 

T 



« 

2.1-4 

2.9-10 

tt 

3 



5.8-9 

1.6-8 

3.1-13 

304 

T 



it 

tl 

4.7-7 

it 

3 
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TABLE II . 10( Continued) 


T e (°K) 

n 

e ~ 

(cnf 3 ) 

e 

Tl 

c 

Line 

(A) 

Case 

# 


IO 1 * 

2.8-4 

4.5-2 

3.4-6 

584 

T 



ft 

ff 

If 

?r 

3 



5.8-5 

1.9-4 

2.8-9 

304 

T 



ft 

tl 

1.9-6 

It 

3 

4 

5 x 10 

io 10 

3.8-8 

2.3-4 

1.6-9 

584 

T 



tl 

2.0-4 

9.0-10 

n 

3 



5-2-9 

1.9-8 

3.4-12 

3o4 

T 



11 

1! 

8.1-6 

ft 

3 


io 14 

3.8-4 

2.7-4 

3.0-8 

304 

T 



It 

5.8-2 

ft 

fl 

3 



5.2-5 

2.7-4 

3.0-8 

304 

T 



It 

IV 

2.7-5 

fl 

3 
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for our assumed Doppler line profile. (For the continuum we require t > Q~ ' 
where £ is the ratio of collisional to radiative recombination.) 

Some representative values of the parameters e and T| (and C) for the 
584 and 304A lines are given in Table IX. 10 . Values were obtained for each 
temperature at n = lO 1 ^ and lO 1 ^ for the optically thin case and case 3 
given in Table II. 3. We can see from the optical depth tables that under many 
conditions the layer is effectively optically thin, i.e. 

1 < t < (e + 'll)' 1 

In this case the total energy emitted in the line from the surface of a layer 
of optical thickness t^ is given approximately by 

E ~ 2 /jt (e +4*-) 1 1 

where (_* = l/B. 

The values of T] and L vary with optical depth and an integration 
over depth is required to obtain the total line intensity from the layer. 

Using the optical depth tables and estimates of 'll and C. from code 1 

we have obtained approximate solutions for some line intensities. The lines 

chosen are the 584, 537, 304& lines and the D3 line. Figures 11.39-3^ give 

results as a function of T for various n while Figures 11.35-39 show the 

e e 

variation with n g for various T g . We note a considerably different variation 

with T and n for the D3, 584 and 3041. line emission. These different 
e e 

variations give some confidence that simultaneous emission measurements 
of the three lines could yield unique temperatures and densities from the 
emitting region. 
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Figure 11-30 TOTAL LINE EMISSION 


n e = lol ° (cm 5 ) 
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Figure 11-32 TOTAL LINE EMISSION 
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Figure 11-3^ 


TOTAL LINE EMISSION 


n = XO 1 ^ (cm" 5 ) 
e v ' 
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Figure 11-35 TOTAL LIME EMISSION 

k 

Te = 10 0 K 









Figure II -56 TOTAL LIKE EMISSION 


Te = 2 x 10 
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4. Population Rates 

It is very instructive to look at the various processes involved 

in populating and depopulating the upper radiating levels. In other 

words, to determine the exact processes contributing significantly to 

the terms TV and 6 in the transport equation. One way to do this is to 

assume a 3 or 4 level atom and evaluate T] and i from the analytic 

th 

solution choosing different levels for the 3 r or 4 levels. This 
would give us an approximate answer. Another way is to use the 
complete solution of the statistical equilibrium equations for the 
populations to evaluate the rates directly into the upper and lower 
levels of the lines in question. 

Solutions of the full SSS equations were used in obtaining the 
relative rate processes shown In Figs. 11-40 and 11-41. Fig. 11-40 
illustrates the processes populating the 584 A line upper level at 
= 40,000°K for various electron densities. Four processes always 
enter significantly. The largest of these is photoexcitation from 
the 2^S level by absorption of photospheric radiation. Direct colli - 
sional excitation from the ground state is next, followed by collisional 
excitation from the 2^S and radiative decay from the 3-4}. Hence, both 
the photospheric radiation and coupling to the triplet levels are impor- 
tant. 


Fig. 11-41 shows relative processes populating the D3 line upper 

level for the same electron temperature and densities. A much stronger 

dependence on n g is noted. Fhotoab sorption of photospheric radiation in 

the D3 line itself is the dominant mechanism at low n g . At high n g , 

this process is small and the collisional rate from the 3"4) is dominant 

3 3 

with collisional excitation from the 3 P and 2 S levels also being 
significant. Thus, the triplet - singlet interaction is again very 
important. 
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Figure II-4l RELATIVE RATES OF PROCESSES POPULATING 
THE d LEVEL OF He I 

Te = 40,000°K 
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III. LIRE TRANSPORT SOLUTION 
A. Basic Equations 

The solution of the line transport equation that we have used in 
the code developed (Code 2) is an iterative solution of the integral 
equation. The details of the method have been developed hy Avrett and 
Loeser (1969) in a manner convenient for the simultaneous solution of both 
line and continuum transport equations. From Section II. D, we note that 
solution for the line and continuum radiation field simultaneously is 
necessary. We will only briefly summarize the method and equations in 
this report. 

In the statistically steady state the rate equation describing the 
pop ula tion of the state i is 


^ (n j p ji - Vi 3 > = = Yji = 0 



- £ 



III.l 


where P. . is the total transition rate from i to j per second per particle 
in the i state. In general, P^ = R^ + C^, where R^ and C„ represent 
the radiative and collisional transition rates respectively. We shall 
assume a Maxwellian distribution for the electrons and helium particles and 
since we also assume a known external radiation field, the transitions 
involving the continuum can be represented by a single term in £q. III.l. We 
can characterize the system of linear equations III.l by a matrix whose co- 
efficients are equal to P^. In representing matrix elements and 

co-factors thereof, we shall always let the first subscript refer to the 
row and the second to the column. Because of the definition of the transi- 
tion rate the subscripts of the P’s will be reversed when substituted for 
the matrix elements a. The general solution of Eq. III.l is 


n. = X P* 1 ; 
1 m 7 


m 


M 

£ ' 
J 


III. 2 
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where is the co-factor of the element a , = P. in the coefficient ma- 
ny jm 

trix represented by the i equations of type mi and N is the total number 
of helium particles per cm^. The matrix of coefficients has the property 
that the co-factor of all the elements in a column are equal, i.e. P 111 ^ is 
independent of m. 


When the medium becomes optically non-thin for certain frequencies 
the radiation field producing internal excitation and ionization for these 
frequencies is no longer merely the external radiation field but is partly 
dependent on the internal properties of the gas and must be determined from 
the radiative -transfer equation 


III. 3 


where cos" (a represents the angle between the direction of propagation and 
the outward normal z and t = fk dz, S is the source function, k is 
the linear absorption coefficient and 1^ is the specific intensity of the 
radiation. In LTE = B^, however in the non-I/TE case must be specified 
in terms of microscopic processes. In terms of such processes the trans- 
fer equation governing the spectral line between upper level u and lower 
level I may be written as 


- 4-jtu -r— = [HtBl, 4hv-nB n ilfhv + 4jtk ] I 
^ dz 1 lu v u ul Y v c \ 


III. 4 


- n A , j hv - 4*k S (T ), 
u ul u v c c' 


where k c represents the continuum absorption coefficient at frequency v q 

and S is the continuum source function. B- , B .. and A are the Einstein 
c lu' ul ul 

transition probabilities for absorption, stimulated emission and spontaneous 
emission. <J v , \|r v and are the normalized emission, stimulated emission 
and absorption coefficients within the line defined such that 


00 CO CO 

S V v = I ♦v dv = 4T J I J v dv = 

O O 4jt O 
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The continuum absorption coefficient is generally very small compared with 
the line absorption coefficient near the line center and will be neglected 
in determining the source function within the line. Using the standard 
relations between the Einstein coefficients and assuming > 

the source function becomes 

„ 2hv 3 1 

s m= — TTiTi^HJTT > 

where g represents the statistical weight. The minus one term in the 
denominator represents stimulated emissions. 


In evaluating the radiative excitation rate R. . for transitions 

cl 

between bound levels the line radiation field enters as 
00 

-f J V ( T ) * V ( T ) dv ■' 

o 

where 

j v ( t) = ET ! I V (T ^ ) 

4 it 

is the mean intensity and du) represents the solid angle. It is thus con- 
venient to formulate the transfer equation in terms of J rather than I . 

v v 

It is now convenient to separate those components involving the unknown 
radiation field, J ul< from the co-factors. This is done by expanding the 
determinant P 1 ^. in terms of its co-factors Q 1 ^. Thus 


P iU = 2 P^ 1 = P-l Q" 1 + 2 P-jjft* 1 

kjfl ^ lu k^l/u 111 


P 111 = z P >Q kU » P n Q lu + 2 P v Q ku . 

k^u ^ ^ k^l ^ 


III. 6 


Actually, may appear in many of the co-factors since the line u-1 may 
fall in the ionization continuum of some other transition. The influence 
of as well as the line radiation in general on the bound-free radiative 
rates will be neglected. Using the standard relationship between the 
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Einstein coefficients, equation III. 6 , and remembering that 

Q 111 = Q lu and J* v dv the source function may be written as 


S Ul P ul pUl 


. f J § dv + eB + L 

piu V V 


1 + e + T\ 


where B is the Planck function 


B v < V " p ul C 


lu 

ul 


e - 


jxl 

\ll 


ul 


2hv 

2 


3 *3 

g. 


u 


L = p , — — T 2 P,v.Q 
^ K ul . _ul , / /, Ik 

A ul Q k ^ u ^ 1 


kl 


T| = 


~"77T S P..^. 


A„,« ul kj^l 


uk 


ul 


III. 7 


III. 8 


The terms entering the numerator of Eq. III. 7 each represent a method of 
populating the upper level from the lower level. The first term repre- 
sents direct radiative excitation, the second direct collisional excitation 
and the third any combination of radiative or collisional processes involv- 
ing one or more intermediate levels in going from the lower to the upper 
level. The denominator, on the other hand, consists of terms indicating 
transition paths from the upper to the lower level. All the terms are 
normalized with respect to A ?1 . The first term represents direct radiative 
de -excitation, the second direct collisional de-excitation and the third 
any indirect process going from the upper to the lower state. 



7^ 

Equation III. 7 is solved using a discrete ordinate method for 
the frequency integral in which we assume 


K 


J F( x )to -s Vx- 

o k=l 


III. 9 


The coefficients are to be found from a given set of dimensionless frequency 
values x^, k = 1,2. . .K; F^ is the value of F at x = x^. The solution 
for is given in Section III.C. 

We shall assume $ to be Gaussian, i.e. 

u 7 


i = 


= exp(-v 2 ) . v _ Av_ 
D 


v " , ^1/2, ’ ' “ Av ‘ 

(«) ' Av.. D 


III. 10 


Avp is the Doppler half width given by 

jo /2kTx lj/2 
c V M ; 

where M is the Helium particle mass. 


III. II 


The source function is obtained at N depth points (i = 1 N) 

within the assumed layer. The depth points are located at specified physical 
depths which do not change during the calculation. The values of are 
obtained from the matrix equation 

N 

£ M S = C. III. 12 

j=l X - J 0 1 

where i and j refer to depth points. 


The coefficients are given by 

_ (A) 


M. 


10 


= A. . - 
ij 


r~T W. . 
1+e . 


where 



I I i 0 = i 
0) j/i 


III. 13 
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W ii^ are functions, e^is a coupling parameter. 


C . are given by 


®i S 1 N r 

C j “ l+e. B i ' IT^r * w ij S j 


s 

In these equations, IL is a coupling parameter. 


and 


ij 


w 

ij 


-jr- E A W. - J ^ k ^1 k 
/n ijk 0 ik +r 


2 K (A) 0-v r - 

= 4“ 2 a w w j fe . j 

/tt . . Tc ijk 0 +r 

J J 


k=l 


where k refers to a specific frequency, A k are weighting functions. 

i, 

profile function 0. is simply 


0 


-V 


ik 


= e 


The optical depth values t are given by 

IK 

T 

T ik = J* ^ik +r i } dT * 


dT 


,L 

k dz 


h.v 

% 3/2 L.% 


n l * 


hi. i4 


III. 15 


III.16 

The 


3 
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The weighting coefficients W . ' represent an expansion of the mean 

1JK 


intensity J^t^) in terms of the source function S v (t ) 


J (t .) = I W.v' ; S (t ) 

v vi' v 

J 


Substituting the solution of the radiative transport equation for 

J we obtain 
v 


1 N /i \ 

tj E ( I t-T I ) S(t) dt = E w \ A) s(T j_) 
11 0=1 ° 


where the frequency subscript has been dropped. To evaluate W. , we assume 

-LJ 

that S . is represented by linear segments between optical depth points n and 
n + 1. S(t) in the interval n £ t ^ T n +q * s £i ven hy 

qi't'i - S ( Tn+1 ) + S ( — ) III. 19 

S^t; - B l - t } n+1 V T - T } 

n+1 n n+1 n 

Equations fcrW.^ are obtained by substituting Eq. IL19 into Eq. II. 18 . 

The ^ so obtained depend only on the set of T. values chosen. Different 
ij 

expressions for are obtained for j > i > j ^ 1 = 1» We give 

id 

these expressions in order. 


W depends only on the coefficient of 
ij n 

t (t-T ,)dt T n+1 

2 w i„= r a =ii^ii t -~T a — + j E ii t - T i> 

t n n-_L t 


t) 


T n+l" T n 


III. 20 
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For n > 1 , 


-It -t . - t -t . 

t.(E 0 |t E_[t -t.|) + e ” - e n 1 

2 1 n-1 i 1 2 1 n 1 1 ' 

~ < 7 n - Vl> 


+EJt -T.I-EJt , -T.l 
3' n i l 3' n-1 1 '- 


[E_ t -t. -E 0 t -t. ] + [E t -t -E t , -t . ] 

t - t , 2' n-1 1 ' 2 1 n 1 1 t ,,-t 2' n l 1 2' n+1 i‘ 

n n-1 n+1 n 


r T i< E 2l T n- T il- E 2l T n+l- r ll >+ e 


- T -T.l T , , -T . 

n I 1 n+1 1 


+ E 3 1 7 n+l” 7 i I ” E 3 1 T n~ T i I 


(l n + l- T n ) 


III. 21 


W" 7 n-1^ r™ /_ _ \ T i^ 


^in “ (r -T n T ) ■ [ E 2 ( T i- T n)- E 2 ( T i- T n-l» + i 
v n n-1 


T , , -T ) 

n+1 n 


[E 2 (T i- T w.l ) - E 2( T i- T n) :l 


+ [e (T± Tn+l) - e ( 1 n) +E 3 (T 1 -T n )-E 3 (T 1 -T n+ 1 )]/(T II+ 1 -T n ) 

- [e ( 1 Tn> - 6 ( 1 n ‘ l) + - E 3 (* !-%)]/( WjT 


III. 22 


For 1 = n, 


2W in = [ 1 - E 2 (T i- T n-X )] + 


[ 1/2 - e ''' i Tn ' 1> + E 3 (T 1 -T n _ 1 )l [ 1/2 - e <Tn+1 ^ 


* V 7 ^- 7 !' 1 


T ~ T 
n n-l 


T n+ 1 ~ T n 


III. 23 
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B. Evaluation of e / and B S 


e‘ and B appearing in Eqs. III. 13 and IIJ.14 are related to the source 
function parameters $ . T| and C as follows: 

e'B S = eB + y, 

111.24 

€ ' = € + Tj 

g 

where B is the Planck function. The values of c 1 and B depend upon the 
levels included in the model other than the upper and lower line levels. 

For two bound levels (1,2) and one continuum level (k) e' and e'B are given 


C P cjd P 

* ' - - e> + x 8i -^ L 

21 21 cu 2 21 


2 21 


where 


2h.^ -hv /kT 

a = — jr > & = e 

c 


P. . = P.. P. . /(P. . + P, .) 

ij lk kj ' ki kj 


For three bound levels and a continuum level 


e ' = X - Y 


B = aY 


III. 27 


where 


(C 01 + P„ n + P^ + P. 


*12^ 


21 21 ‘23 23 Mgg 


X 


III. 28 
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Y 


21 


- (C 12 + P L> 


M n „ R~ 

12 2 


*22 


) 


*12 = P 32 + P 32 

= P 23 + P 23 

*22 = P 31 + P 31 + P 32 + P 32 

R 2 = P 13 + P 13 


III.29 


XII . 30 


For some applications it is convenient to solve for two line radiation 
fields simultaneously using only the four levels in Code 2. The two lines 
are then represented by the 3 ~ 1 and 2-1 transitions. The corresponding 
X and Y values for the 3 ~ 1 lihes are 


l - - M 1P 

X - AT < C 51 + P 51 * P 52 + p 32 - S^ 1 ? 




0) 

X ^ + ^3 “11 


where 


% 


p + P 4- P 4- P 

21 21 r 23 r 23 


R 1 P 12 + P 12 ' 


III. 31 


III. 32 


In the general ease e ' and B S are obtained from Code 1 at each 
optical depth. Code 1 can be run for arbitrary values of J for the lines 
or A for the continuum rates. 



8o 


c. 


which 


where 


when 


when 


where 


Evaluation of 

The evaluation of A^. follows that of Avrett and Loeser (1969) 
is recommended for further details. F^,(x) represented by 


K 




0 < x <■ 


*K 


f(x) . = 1, 0 £ x £ x^. 


j = I, and 


f(x) j = 


(1 - — ) (1 - y — )> 0 ^ x ^ x 

x o A 5 


0 , Xj £ X £ x^ 


) = 2,3,...K.ly is an adjustable parameter. 


III. 33 


III. 34 


is given by 


K 




= 2 g, f 
j=l 


-1 


j jk 


is the inverse of f(xj 5 .)^ and 


g j = 




*K ' ^ ‘ 1 

2 ^j ^ = . • .K. 


III. 35 


III. 36 
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D. Sample Solution - Two lines 


The Inputs to code 2 to begin the iterative solution are the 

n^/n^. Approximate starting values are obtained from code 1. We shall 

illustrate the solution by discussing the populations n^/n^, n^/n^ and 

n^/n^, where 1, 2 and 3 designate the lower level and two upper line 

levels for which the transport equations are solved and K the continuum 
s 

level, and are obtained as described in Section III.B. n^ and 

the optical depths at the geometric depth points are calculated by the 

code from the given population ratios and For the initial mm, 

the upward radiative rates are those given by the optically thin case. 

Sj values are calculated from Eq. III. 12. is found from the S^ and 

the new radiative excitation rate based on J . is used to obtain new 

J 

values of n^/n^ through code 1 or a simpler 3 or 4- level solution. The 
new values of n^/n^ and ** are used as input to code 2 (or 3) to obtain 
the other line intensity. 


We illustrate the solution for two lines in Fig. III-l for T g = 

2 x 10^ °K and n g = 10 12 cm"^. The Figure shows results for n 2 /n , 
n 3^ n l an< * n K^ ri l aS a ^ unc ^^- on ^he number of iterations. The initial 
values of n^/n^ and n^/n^ are 10 - ^ and r^/r^ is 10 2 . The solutions 
have essentially converged after 2-3 iterations. Solutions at the surface 
of the layer are labeled -0; solutions at the center of the layer are 
labeled -C. The codes have been tested for a number of simultaneous line 
and continuum transfer problems. There have been no convergence difficul- 
ties provided the lines and continuua chosen for solution do significantly 
effect each other. 
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IV. CONTINUUM TRANSPORT SOLUTION 

A. Basic Equations 

The statistical equilibrium equation for one bound level ( n^) and 
one contiuum level (n. ) is 


s [ f n - 

1 v x 


"k 


CO 


-hv/hT 


wv ** y OVt 

^ [ ^ I 6 ( + J v ) dv + C ] , 

Tc v n c 


IV. 1 


where the * refers to the LTE population at a given electron temperature and 
density. It is convenient to express the population ratio n^/n^ in terms of 
the non-equilibrium parameter b^. 


b l = / (“k/ 13 ^*) * 


IV. 2 


The radiative transport equation for v > can be written 



IV. 3 


where neglecting stimulated recombination (which is Justified for He for the 
temperatures considered) 


B 



-hv/kT 

e 


dT 


v 


k dz 
v 


k = n- a 
v 1 v 


IV. 4 
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The values of b^ at depth i are found by solving the set of equations 


N 1 

Z M. . -r — 

>1 b ij 


e i 


e i 


i = 1, 2 N 


IV. 5 


where 


M 


i<5 




e? R» 

1 X 


A. . = { lf J " 1 

0, 3 + i 


K 

2 

k=l 


A’, W 


ijk 


( A -D fk' e ’ y k V 


IV. 6 


and 


R*. = 


K 

2 

k=l 


Slj 


\ r 


■V 6 i 


IV. 7 


is the optical depth at the continuum threshold at depth 1,0^^= h.v^/hT^, 

= v^/v^, k = 1, 2, . ,.K, a set of dimensionless frequency values and 
are values such that g^. y k ^ represents the frequency dependence of the photo- 
ionization cross section. The optical depth is given by 


fk 1 

«i J7 - r 


IV. 8 


s. "b 

and are coupling parameters whose calculation is described in the next 
section. The values of are related to as follows. 


„ (A-1) _ , V ifr (A) - X ’ 3 * 1 

i3k v. „. ( A) , J / i. 


IV. 9 


ijk 


The mean intensity of the radiation field is given by 


j g w ( A) B Jk 

« j=l 13k b ij ’ 


IV. 10 
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Evaluation of e 


£L Id 

€ ^ and follow from the statistical, equilibrium equation solution 


for b^. 


, . a 
1 + e 


IV. 11 


R-. = 4* f —• J'dv 

Ik J hv v 


IV. 12 


co a 

v 


Si* -M 57 ^ 


For two bound levels (1,2) and a continuum level k, e^ 3, and are 


e b - — 


) 

* lk V 2k 

Ik 


a 1 , r ^ 2 /“l ) P 21 

5 = * (C lk + P +P 

R, , 21 2k 

lk 


* , *. 


21 ^k / n 2 ) P k2 


IV. 13 


For an N level model the corresponding equations are 


b 1 ^ 1 1 

6 “ R '^lk kl ^ . i 

R ik lk (T- l ± k 


p u 


s a --i-(0 -f-i- E 

®Ek 111 a 1 * I 


Q A fl J* k 


P 

Q / ■ 


IV. 14 


a b 

€ and e can easily be obtained from program 1. 
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C. Evaluation of A / 

k 


The 




coefficients are defined such that 


y i 

J* F(y)dy = S A, F IV. 15 

J 1 k=l 15 k 

where F, is the value of F at y = y . Here, 1 < y £ y . With unity 
K K K K t 

instead of zero as the lower integration limit, we determine the A^ 
coefficients as follows. 


Let f’ (y)^ =1, 1 ^ y ^ y k . 


when j = 1, and 


f (y)j = 


(i- J4i> (i - * fri). i* 

y j y j 


r 

' 0 , 


IV. 16 


y , ^ y ^ yi 


J 


K' 


when j = 2,3>...K. Let f^”^ he the inverse of f'(y k )j. 


Then 


K 


^ "j?! S ' J f ' 


IV. 17 


where 


ri - 

3 j 


{ 

■ o 


y K -i> j = i 


iv. 18 


| (y 3 -l)(l- |), j = 2,3,, 


.N. 


Y is an adjustable parameter. 
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V. CODES 

A. Code, JL -_Solutlpn of .Statistic ally Stea dy State Population 

Equations 

1. Main Program 

In our own internal notation, this code is designated as P4-8. The 
main program reads f numbers, wavelengths (lambda), electron and radiation 
temperatures, the energy level model, recombination coefficients, con- 
verts energy differences from wave numbers to electron volts, and writes 
all the above before entering subroutine CONSTS. The subroutine CONSTS 
as well as other subroutines are described in the next section. The 
radiation temperature is that characterizing the external Blackbody 
radiation field incident upon the layer. The indices on the f numbers 
refer to the ion and the upper and lower line levels respectively. 

On return to the main program, more quantities are read in and 
printed out, such as: the total line center optical depth of the layer, 

the dilution factor for the external radiation, the electron density, 
the indeces of the chosen thick line, the depth points at which 
solutions are to he obtained, a constant divisor for the elements of 
the main matrix to prevent overflow, Y values for thick lines, and 
various control numbers. These quantities are entered here in order 
to be able to run a series of different solutions by changing any or 
all of them without changing the electron or radiation temperatures 
and, thereby, having to recalculate the many reaction rates. 

After the above have been printed and the ratio of Bq/E^ cal- 
culated the subroutine AMAT is entered to calculate the elements^nd 
cofactors of the A matrix for the optically thin case. The rate equa- 
tion matrix represented by Eq. II- 3 is printed out. The matrix is 
30 x 30 corresponding to the number of energy levels. Selected cofactors 
are printed out. The accuracy of the solution can be checked by compari- 
son of cofactors with the same j values. 
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The solution for the non-equilibrium populations n. is obtained 
in AMAT. The solutions axe printed for the optically thin case along 
with LTE solutions for the same temperature and density. The optical 
depth at the center of each line is calculated and printed out for a 
layer of optical thickness T 1 in the line specified in the input. 

Some further quantities such as the physical thickness of the layer H 
are calculated and the subroutine ELIM is entered to compute e, T[ and C 

The program now repeats the calculation with the specified set 
of Y values read in. Next, a IX) loop on depth is set up. In the list- 
ing given there is no depth variable. In calling for one to obtain 
values of T) and C for example, it would enter into the calculation of 
the quantity FF, which Is used to alter the Y value for the thick line. 
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2. Description of Sub -Programs 

Subroutine AMAT(KK) 

AMAT calculates elements of the A matrix using quantities determined 

by CONSTS. Once all the elements are found, they are divided by a constant 

to prevent overflow. Cofactors are calculated, summed, and used to determine 

X. From this the n. are found. Next the n . are found and the ratio n./n .. 

i eqi l eqi 

Last, before returning to the main program, the optical depths are calculated 
and printed for each f value. 

KK = 0 , optically thin solution 

KK = 1 , optically thick solution for any Y 0.0. 

Function COFACT (NR. NC. ME. D) 

COFACT finds the cofactor of matrix D for row NR, column NC. NE is 
the number of rows and columns in D. 


Function CflLL (T. I. J. K. N> 


Function C0LL calculates certain electron collisional excitation and 
ionization coefficients for special points. 


T 

I, J, K 
N 
N 


electron temperature (°K) 
indeces of coefficient 

1 , excitation coefficient calculated 

2 , ionization coefficient calculated 


Subroutine CONSTS 


Initially, if the radiation temperature is different from the electron 
temperature, CONSTS reads in new values of p_ k calculated at T r . It then 
calculates the electron collisional ionization and recombination coefficients, 
the radiative ionization rate, the B„ k and Et^., the electron collisional 
excitation coefficient, the Einstein spontaneous transition probability, and 
l.ho Einstein absorption transition probability multiplied by the Planck function. 
The subscript i refers to the ion, the second subscript the beginning level . 
and the last subscript the ending level of the transition. 


2 
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AP = Einstein spontaneous transition probability for the thick 
line (K7, K8, K9) 

RHO ss hv^/c^, where v is calculated at (K7, K8, K9) 


This computes double cofactors for use in calculating eta and iota. 
It is called from the main program and uses the function ELM2. 


Function ELM2 (IRA. ICA. KRA. KGA. RMAT. 


IRA = row index of 1st row to be eliminated 

ICA = column index of 1st column to be eliminated 

KRA = row index of second row to be eliminated 

KCA = column index of second column to be eliminated 

RMAT = contains reduced matrix 
Z - sign of cofactor (Z = ± l). 


ELM2 eliminates two rows and two columns from the A matrix and places 
the reduced matrix in RMAT. 


3 
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B. Radiative Transfer Codes 

Two codes have been developed under this program. The first 
(code 2) represents a numerical solution of the line transport equation 
given in Section III. The second (code 3) represents a numerical solu- 
tion of the continuum transport equation. Listings of these codes appear 
in Appendices B and C. The simultaneous solution of several line and/or 
continuum transport eqs. with the population eqs, is accomplished by- 
iteration. Initial values for population ratios are estimated from 
various solutions obtained from code 1 (see Section II). After obtain- 
ing the radiation intensities from codes 2 and 3 the new radiative 
rates are used as input to code 1 to obtain new population ratios to be 
used in the next iteration. 

Code 2 requires the specification of certain quantities involving 

t g o. "b 

e.^ and while code 3 requires and e^. There are two methods provided 
for determining these quantities. One option is to calculate them in 
code 1 and simply read them in. (Code 1 is now set up to calculate the 
quantities 1] and (, which are read in and are related in a simple manner — 

Q 

see Eq. III. 24. The quantities e and e are easily related to the co- 
factors generated by code 1. Code 1 will generate arbitrary cofactors 
depending upon input parameters. ) 

A more approximate method of calculating these quantities is pro- 
vided directly in the two codes. Code 2 will obtain the solution for 

e , and B from a 4 level model (3 bound levels and 1 continuum level). 

9 > t 

Code 3 has the option of solving for e and e from a 3 level model 
(2 bound levels and 1 continuum level). These approximate solutions can 

1 

be very useful for some purposes. 

A symmetrical (about the center of the layer) geometrical grid 
(called Z in the program) is set up given the total geometric depth and 
number of desired decades. Currently, the maximum number of decades is 
set at 5 giving a total number of points of 31. This calls for solution 
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of a 3l x 31 matrix and may create underflow or overflow problems. 

The number of decades required for convergence must be tested for each 
problem. 


Optical depth points at each geometric depth are next calculated 
from an integration over geometric depth at a given frequency. 


x = 


v- v 


Av r 


Steps in frequency are now calculated using the quantities 
in the bound level cases and y = vf in the continuum. 


Step 


sizes are input quantities and are constant in x and y. The maximum x 

—x^ 

is found by choosing the first value of x such that e t ^1.0, 

^ max 

and the maximum y such that t /y < 1.0. 

max 


Certain quantities uniquely dependent on frequency are now calcu- 
lated, such as f(y) , g(y), A(y), f’ (x), g' (x), and A' (x). A loop over 
frequency is set up with index k to compute the W and W. ^ ^ (by 

ijK m jK 

means of the subroutine WMAT) and then sums over frequency are taken within 
the loop to arrive at the coefficients in the final equations: W„ in 

the hound case and M. . in the continuum. 

Once the integration over frequency is completed, the program 
sets up a matrix, called EM, for solution of the simultaneous equations 
to get S. or l/h , . . These are then punched in cards along with the pop- 
ulation ratios to be used as inputs to the next step in the iteration. 


The S , or b . . are used to determine new line or continuum radiative 
<5 ij 


rates which are then used as input to code 1. New values for the popula- 
tions and the parameters B®, and are obtained and are used in 
beginning the next step of the iteration. 
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DESCRIPTION OF SUBPROGRAMS 


SUBROUTINE NMAT (HZ, TOUT) 

WMAT calculates W. ^ and W for those areas of the matrices 

■LJ ij 

where is less than TCUT. 

FUNCTION NOSONI (A, X t L, IldAX) 


NOSONI is a matrix inversion routine using the method outlined on page 
434- of Hildebrand, Introduction to Numerical Analysis (New York, 1956). 

A is a matrix of order L with column dimension LMAX. Its elements are 
assumed to be stored columnwise in the usual Fortran manner. X is work- 
ing storage of length L. The inverse of A will replace A. Upon return, 
NOSONI = 1 if inversion went properly, = 0 if a divisor is zero, in which 
case, A may contain garbage. 


FUNCTION ESB (X, N) 

This routine uses the function EXIN to obtain E^(x), the exponential integral 
of order 1 in X. It then calculates Ejy(X), where N = 2, 3 or 4 from the rec 
recurrence relation E^ + ^(X) = i [e X - ZE^(X)] 

FUNCTION EXEN (Y) 

EXIN obtains E^(Y) from polynomial approximations shown in sections 5*1*53 
and 5*1*56 of the Handbook of Mathematical Functions , National Bureau of 
Standards, 1964. 


I 



94 


VI . REFERENCES 

Allen, C. W. : 1961, Mem. Soc. R. Sc i. Liege, 5th Sec., 4, 24l 

Allen, C. W. : 196?, "Astrophysical Quantities," 2nd Ed., The Athlone Press, 
London, p 4l 

Athay, R. G. , and Johnson, H. R.: i960, Ap. J. , 151. 413 

Avrett, E. H. and Loeser, R. : 1969, SAO Special Report No. 303, Cambridge, Mass. 
Benson, R. S. and Kulander, J. L. : 1972, Solar Physics, 2J_, 305 
Burgess, A. and Seaton, M. J. : I960, M. N., 120, 9 

DeJager, C. and DeGroot, C.: 1957; B.A.N., 14, 21 

Dolder, K. T., Harrison, M. F. A. and Thonemann, P. C. : 1961, Proc. R. Soc., 

A264, 367 

Green, A. E. S.: 1966, AIAA J. h, 769 

Hearn, A. G. : 1969, M. N. , 1 k2_, 53 

Jefferies, J. T. : 1955* Aust. J. Phys., 8, 335 

Jefferies, J. T. : 1957; M. N., Ill, 493 

Kieffer, L. J. and Dunn, G. H. : 1966, Rev. Mod. Phys., ^8, 1 

Long, D. R. : 1967; Thesis, Univ. of Wash., Seattle, Wash. 

Miiialas, D. and Stone, M. E.: 1968, Ap. J., 151, 293 

Peach, G. : 1967; Mem. Roy. Ast. Soc., 21, 15 

Saraph, H. E. : 1964, Proc. Phys. Soc., 83, 763 

Seaton, M. J. : 1962, "Atomic and Molecular Processes" (ed. by D. R. Bates), 

Academic Press, N.Y., p 4l4 

Shklovsky, I. S. and Kononovitch, E. W. : 1958; Russ. Ast. J., 57 

Stewart, A. L. and Webb, T. G. : 1963, Proc. Phys. Soc., 82, 532 

White, R. 0.: 1961, Ap. J., 1^4, 85 

Wiese, W. L., Smith, M. W. and Glennon, B. M. : 1966, Atomic Transition 

Probabilities, NSR DS-NBS4 

Zirin, H. : 1956, Ap. J. , 1 2J, 536 



APPENDIX A 

PROGRAM FOR SOLUTION OF STEADY STATE 
EQUATIONS 


! * 
2 * 
3* 
a* 

5* 

6 * 

7 * 

e* 

9* 

JO* 

M* 

1 2 * 
13 * 

J 4 * 

15 * 

1 ** 

jT* 

1 «* 
1 9 * 
20 * 
21 * 
? 2 * 
?S* 
?4* 
?5* 
2 ** 
27 * 
? 8 * 
29 * 

30* 

31* 

32* 

33* 


C P • 48A HEUUM 

c single precision 

c VARIABLE r«*zero 

C ALL BETA AT TR 

DIMENSION DA(2) 

OyMENSlON EV<3«19*19) 

,CoMMON/8LKt/TAU(iO) »XJ(10»4) »NTAU 

J /BLK2/0MEGAC3* 19» 19) t aP ( 3» 19» 19) »0B(3» 19* 19) ,nBAR(i. 19. lo> 

(j/BLKfc/LEfll.J^iafK^KRUROiICO 

J^2» (4) » Tl*$lBETA*SlCfSlAPtSlR7 

8/BLK8/zAN(40),XNt0(4 0 ) 

9/BLKl3/NL»NP»NT»MA(3) * N5T ( 3) *MI ( 10 o) * MJ ( 1 QO) . mK M O n) » NALPH 

C O MMON/8LKl5/NrlCl00),NFj(100) f NFK(iJ 0 ),NOF '** 11 * { > UH * LP * 

CALL DaTE(9»0A) 

I?YX « 15 
10 DO 15 1*1*00 
DO 15 J*t*«0 
15 A(I*J)sO a 
DO 20 1*1*3 
DO 20 J* 1 ♦ 19 
DO 20 K*1 t 19 
BfTACIf J»K)«0# 

EZ C I* J*K) a 0,0 
ABARU, J*K)«0, 

Y(I*J*K) * 0.0 
AlPH(I*J.K)«0 # 

AP(IiJ,K)*0. 

BB(I*J f K)«0. 

0BARCI,J*K)*0, 
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3«* 

35* 

56* 

37* 

38* 

39* 

40* 

at* 

42* 

43* 

44* 

45* 

at* 

47* 

48* 

49* 

so* 

51* 

52* 

53* 

54* 

55* 

5t* 

57* 
58* 
«jt* 
60 * 
61 * 
62* 
63* 
64* 
65* 
6t* 
67* 
68 * 
69* 
70* 
7 1 * 
72* 
73* 
74* 
75* 
76* 
77* 
78* 
79* 
80* 
81* 
82 * 
83* 
84* 
85* 
86 * 
87* 
88 * 
89* 
40* 
61 * 


D0AR(I* J*K)*0, 

OMEGA(l»JtK)*0,0 
20 F(I*J,K)*0. 

READ (5*40) NF, NT »NALpH*NOF 
Nl s NF t NT t l 
25 F0«MAT(3I6,E12,8) 

00 30 !I*1*N0f 

read (5*?5)I*J,K,(P(j,J»K)) 

Hl(II) « I 
MJ(II) * J 
MK ( 1 1 ) * K 

nfKII) • I 
nf J (II) * J 
NFK(II) 3 K 

30 continue 

35 FoRMAT(3jt f 15 x*2E15,8) 

40 F0RMAt(12!6) 

45 FORMAT(6Cl2»0) 

READ C5*45)T*TR 
XlAM s' 5,876 E»o5 
ANU = Cl/XIAM 
IF (TR»1.0) 55,55*50 
50 EX * (H*ANU)/(XK*TR) 

DEN * EXP(EX) - 1,0 
By » £2.0 *H*ANu**3/Cl** 2)/DEN 
Go TO 60 
55 BV » 0.0 
60 CONylNuE 

65 FoRMAT( 28H1ELECTRQN TEMPERATURE <TE) MPEtl , 4 , lOX*j?8HRAOlATlON 1 
IPERATUrE (TR) sE 1 1 .4, 10X , 5HBNU *E1 1 .4) 

Do 70 1*1*3 
DO 70 J* 1 * 1 9 
D 0 70 K*l*l9 
IF fJ.LE.K) GO TO 70 
F(I*K,J)*F(I>J,K) 

70 CONTINUE 

75 FORHAT(7HlHlLlUM*10X,4HP48A*lQXf8HPART ONE* 70 X.a6 # A3) 

WrITE (6, 75)0 A 

80 FORHAT(9HOF VALUES, 2 a X*6HLAH8DA/) 

WRITE (6,80) 

85 FqRNAT(SH F(I2,IH,I2,iH,I2,3H) = IRE 1 0 * 3 ,E15 . 4 ) 

DO 90 1 1*1 , NQF 

1 * HI ( I I ) 

J * MJ(II) 

K * MK(II) 

XL * 1 . 0 E* 08 /A 88 (E(I, J) • E(X*K)) 

WRITE (6,85)I*j*KtF<!.J*K),XL 
90 Continue 
00 95 J*1»100 
MlCJ) * 0 
Mj(J) a 0 
MK ( J ) a 0 

95 continue 

do 100 ii*i*nalp h 

READ (5*35)I*J,K*(ALpH(I*J,K)),8ETA(I*J*K) 

Ml(II) * I 
MJ(II) * J 
MK(II) * K 
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92* 
93* 
99* 
95* 
99* 
97* 
98* 
99* 
100* 
lpt* 
1 02 * 
10 3 * 

1 04 * 
105 * 
lo8* 
107* 
toe* 
io 9 * 
MQ* 
ill* 
M2* 
M3* 
M4* 
M3* 
I j6* 

M 7 * 

M»* 

ll9* 

120 * 

121 * 

1?2* 

1 2 3 * 

124* 

123* 

1 |8* 


l?8* f 

129* 

130* 

131* 

132* 

133* 

13«* 

135* 

13 ®* 

13 7 * 

138* 

139 * 

140* 

iai* 

142* 

143* 

1 il4* 
145* 
148* 

U7* 

I48* 

1«9* 


100 CONTINUE 

W&ITE (6,65)T,tR»BV 
WRITE (8,105) 

105 FO«HAT(1HO,22X,11HSTaTI 8TICAI,15X*9HNUMBER Op, 17 Xm 9h ENErGy VALUES 

» aB0VE/3X»3HI0N*3x,5hLEVEL»11X,6H*EIGHT»12X*21H0uTfR SHELL ELECTRO 
?N8t9X,22HE(lM) IN W*VE NUMBERS) 

DO 110 !sl.3 
NJ R NsTCl) 

WRITE C8,H5)l,CJ,G(l,J),ZBCI,J),E(I,J)tJ«l,NJ) 

110 CONTINUE 

115 Format ( 1H0 * 14, l7»Fl7.0»E25,0»F34,2/(I12fP 17 , o,F25,o*E34*2)) 

120 FQRMAT(iH0.1ZX,3Hl », 2 X l iHg. 1 9X , 1«3//4X »lMK*gX* 3HJ «»?X • |H| » 19X t lM 
11 /) 

125 FORMAT ( 15 » lP2E20«7) 

130 FoRMAT( 8 H 0 l * 1 ) 

135 F0*MAT(lH0,3X, 1HJ*8X,3HK «.2X» 1HM 17X* 1H2* 17 X , 1 h 3.|7X, 1H4» l7X , 1H5. 
1 17XtlH6»l7XtlH7/) 

140 FORMAT (1M0,10X,3RK ■,2X,lHS,l7 X ,lH9,i7X,2Hio #1 6X,2Hll,l8X,2Hl2,18X 
i#2H|3,l8X,2Hl4/) 

145 FORMAT (1 MO, 3X,1HJ»8X,3HK « , 2X * 2 H 1 5 , 1 6 X ,2H1 6 , 1 &X ,?Hi 7, 1 6 X , 2h18, 18 X i 
12H19/) 

150 F 0 RMAT( 8 HQI * 2) 

155 FORMAT ( tM0,3X, 1 HJ, 8 X, 3HK »,2Xt 1 HI ♦ 19 X » 1M2. 19*, 1H3. |9X, IH4M*X, 1H5, 
?/) 

180 FORMAT (1 HO, Ux, 3HK .*,.2X*|H6f I9x« lM7f 19X*1M8,19X, 1M9,19X*2H10/) 

165 F0RMAT(l5»lP7Et8,7) 

170 FORMAT (I5«IF5ej8*7) 

175 F0RMAT(I5» 1P5E20.7) 

180 F0RMAT(iSHlALPHAM»J*K)*26X*23HREC0MBlNATl0N COEFFICIENT) 

WRITE (6,180) 

WRITE (8,120) 

WRITE (6,l25)(Kt*LPH(2,l,K),ALPH(J,l,K),K*l, 1 9) 

DO 185 1*1*3 
JL * NST(I) 

DO U5 J*l,JL 
DO 185 K*t » JL 

E 2 CMJ,K) * A 0 8 (E(i,j) . 'E(ItK)) 

EV(X*J,K) * i.23977E»04*EZ(I»J,K) 

185 CONTINUE 

WRITE ( 6 . 190) 

190 P0RMAT(t0HtDE(i»J,K),29X,36HENERGY DIFFERENCES IN FLEcTRqN VOlTS) 
WRITE (6,130) 

WRITE (6,135) 

WRITE (6,1 65) ( J , CEv ( i , J t K> , K*t ,7 ) t J* 1 ♦ 19 ) 

WRITE (6,140) 

WRITE (6,16S)(Jf (EV(t,JtK),Kii8,l4) ,J*1.19) 

WRITE (6.190) 

WRITE (6,130) 

WRITE (6,145) 

WRITE (6,t70){J,(lV(i,J,K),K«iS,l9),J»l,l9) 

WRITE (6,150) 

WRITE (6,155) 

WRITE (6.175) (J.(EV( 2 . J»K) ,K«1 ,5) * J*M 10) 

WRITE (6,160) 

WRITE (6,175)(j,(EV(2,J»K),K*6,i0),Jpl»10) 
c calculation Of constants 
call consts 

195 CONTINUE 
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150 * 

151 * 

152 * 

153* 

154 * 

155 * 

I56* 

157 * 

158 * 

15 ** 

1 * 0 * 

161 * 

ua* 

1*3* 

164 * 

165 * 

166 * 

167 * 

1 * 8 * 

1 * 9 * 

170 * 

1 7 1 * 

172* 

173 * 

174 * 

175 * 

176 * 

177 * 

178 * 

179 * 

lao* 

181 * 

laa* 

103 * 
184 * 
185 * 
106 * 
187 * 
188 * 
189 * 
190 * 
191 * 
192 * 
193 * 
194 * 
195 * 
196 * 
197 * 
198 * 
199 * 
200* 
?0l* 
20?* 
2o3* 
2o«* 
2ft5* 
206* 
26 7* 


READ (5»45)T1»w»XNE 
READ (5»40)K7 f K8*K9fL.AST t UE»NTAU 
READ (5»45)(TAU(t)»I«l»NTAUJ 
READ (5* 45) DIVIDE 
200 READ (5#40 )NOy 

READ C5*40)IPROB*IUsT 
Do 205 1 * 1 *3 
00 205 J * 1 • 1 9 
DO 205 K* i 1 19 
V C I » J»K) « W 
205 CONTINUE 
NSWT » 0 
WRITE (6,2103 

210 EoRNATf lNl*53x f 7Hp . 40A/54 X, HHsInGLE PrECI S ! 0N/54X » i'S'HyARl ABLE 
i*ZER0/54X#27HGAUSS • HERMITE QUADRATURE) 

WRITE (6*2153 TltW*XNE#K7*K8»K9,UE»NTAU»0lVlD E# N8WT*(TAU(I)fI«ltN 

Tuj 

215 Eo*MAT(5H0t 1 *1PE10.3/4H W *El0.3/6H XnE *E1o.3/5H k7 *I3/5H K 8 
13/5H K9 * 1 3/SH It *Is/7H NTAU *I3/9H DIVIDE a£10,3/7H NSwT *13/6 
?TAU «(E10.1)) 

WRITE (6,220)NF,NT f Nt 
220 FORMAT (5H NF *I2/5H nT *I2/5H NL *12) 
write (6,225) IpR&B 
225 F0PRAT( ISHl PROBLEM NUMBER* 143 
If (NOY) 230 *260*230 
230 WRITE (6,235) 

235 FQRMATC9H0Y VALUES/) 

240 F 0 RmAT(1H Y(I2,1H,I2,iM*I2,3H) -1PE9.2) 

DO 255 U*l*N0y 

read (5*25)i*j,k*cy<t*j*k3) 

XNU * CL*Ez(I,J,«) 
ey * (H*XNU)/(XK*T) 

EyZ * (H*XNU)/(X r *TR) 

if ((E¥,eT,e5, 0 )»OR,(EYZ,GT*85.0)) GO TO 245 

BO ■ 2.0*H*XNU*<XNU/oL**2)*XNU/(EXP(H*XNU/{XK*T 3) - 1.0) 

8 T R » 2,0*H*XNu*CXNU/CL**2)*XNu/(ExP(H*XNU/(xK*TR3) *• 1.03 
B08 TR * 90/BTR 
GO TO 250 
245 continue 

exY ■ XNU*(H/XK3*(1, 0 /TR * t.O/T) 

80BTR * expcexy) 

25o Continue 

Y (I * J*k3 » Y(X,J»K)*ro8TR 

YU*K*J3 * Y( 1 , J*K) 

WRITE (6*2403I,J*K*Y(I*J*K) 

255 CONTINUE 
GO TO 275 
260 Write ( 6 , 2653 
265 FORMAT ( 17HqALl Y(I*J,K> 1 W) 

WRITE (6.2703 

270 FORMAT (25H ALL A0AR MULTIPLIED BY W/37M aLL BB<lW,K3 MULTIPLIED 
tY Y(I,J*K)/30H B ( I f J,K) IN BB TABLE IS AT TR) 

27S WRITE (6 f 280)XNE,W 

280 FORMATfSHtNE * tP£ H * 4 , 1 OX * 3HW sOPF7,23 
WRITE (6*205)K7*X8*K9 
285 FORMAT (1H03X2HK7*4X2hK8,4X*HK9/3I6) 

CALL AMAT(O) 

HEM a 6,69E«24 
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20** 
2 q 9 * 
?1 0 * 
2 jl* 

2i 2* 

213# 
2l9* 
2l5* 
2l 6* 
21 7* 
2 1 8 * 
2 1 9 * 
2 ? 0 * 
221* 
2?2* 
2?3* 

2 ?a* 

225* 

2?** 

2?7* 

2 ? 8 * 

229 * 

230* 

231* 

232* 

235* 

23«* 

235* 

236* 

237* 

238* 

239* 

2a 0 * 

2fll* 

2a2* 

2a3* 

2a4* 

2a5* 

?« 6 * 

2<l7* 

2<|8* 

2i|9* 

2 ? 0 * 

2S1* 

253* 

253* 

25 «* 
255* 
256* 
257* 
258* 
259* 
2ft 0* 
2ft 1 * 


SQ 8 2,0*XK*T/hEM 
XnUM * SORT(SQ) 

Dv * XnUH**NU/CL 

A9U * 7*o6E*07 

SlG « 3.79E*U/SQRT(t) 

17 * MACK7) 

0* b T 1/ CZAN C 173 *SIG> 

Of * ZAN(9)*A9fl*tH*ANU)/DV 
Of ■ (PX*0X>/(4»0*P!«SQRTCPX)) 

WrItE ( 6 « 29 0 ) D I * 0 V » Dy 

290 ro«MAT(9H0I/OELT «1PE1 1 »4* 1 0X»aHDV sE 11 • 9 t 1 0* , UMflX *Ell,a) 

KK * 0 

write (6,295)iro»ico 

295 FoRMAT( I 9 H 0 EEIMINATE ROW, 1 3 » 5X , 6 HC 0 LUMN , 13) 

XnU»CL*EZ<K7*K8*K9) 

RHO * (XNU/CL)**2*2,o*H*XNu*GCK7»K9)/G(K7tK8) 

EP5*0MECA(K7tK8»K9)/AP(X7*K8,K9)*XNE 

El«E(K7«K8j 

E2=E(K7»K9) 

CALL ELlMtRHO, aP(K7»k8*K9J) 

0 0 325 LV*l»NTAU 

EY * CM*XNU)/CXK*T) 

F.yZ * (H*XNU)/(XK*TR) 

IF (CEY.GT*85,o)iOR, ( EYZ,GT.e5,0)J GO TO 305 
B0 B 2,*H*XNU*<XNU/Ct**2)*XNU/(EXP(H*XNU/(XK*T) W.) 
H T R*2 .*H*XnU*(XNU/CL**2)*XNU/(EXpIH*XNU/(XK* tr )),1.0) 

RqBTP s 80/BtR 
GO TO 310 
305 Continue 

ExY XNU*(H/XK)*(1,0/TR * 1.0/T) 

0O0TP 3 EXP(Exy) 

310 CONTINUE 
FF * BpBTR 
Y(K7»K8»X9) * FF 
Y(K7fK9fX8) * FF 
B0C*BBfK7*K8»K9)*FF 
BrD*0B(K7*K9»K8)*FF 

WRITE (6«3l5)TAGiLY),K7»K8,K9»BBC»K7»K9»Ke»BRDfK7»K8»K9*K7*K9tK8» 
IF 

315 FoRMATfftHtTAU c 1 PE 1 3 .6 * 9 X , 5HBB* Y ( 12 , i H, 12 , \ H, 12, 3 H) =El3,6,«X,5HE 
1 *Y(l2*lH»I2tlH*l2*3H) *t 1 3 , 6 * <lx f 2HV C l2 ♦ 1 H* 12 , t H* l2 , 9h) , 

?lHtI2»3H) SE13.6) 

WRITE (6 f 320) 

320 F 0 RmAT(89X,13Hy Is CONSTANT) 

CALL AMATU) 

OX * T1/<*AN(I7)*SIG) 

01 a ZAN( 9 )*A 9 i,*(H*ANU)/DV 

Of * (DI*1>X}/U,0*PI*SQRT(PI)) 

WRITE {6t290)Dl»0v»Dy 
CALL ELlM(RH0tAP(K7*K8*K9)) 

325 CONTINUE 
33ft CONTINUE 

GO 70 (10tl95»200),I L AST 

EnD 


)lARNOSTfr S 

ktfON TIMF s «, 79 CPU SfcCONUS 
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1* 

2 * 

3* 

a* 

5* 

b* 

7* 

8 * 

9 * 

10 * 
1 1* 
12* 
13* 

|4* 

|5* 

1 7* 
1»* 
1 9 * 
? 0 * 
?1* 
22* 
?3* 
?4* 
2** 

2** 

?7* 

20 * 

? 9 * 

30 * 

31* 

32* 

33* 

3** 

35* 

38* 

37* 

30* 

3’* 

qO* 

at* 

« 2 * 

a3* 


SUBROUTINE AMat(KK) 

dimension cap ( 40 * 40 ) ,riuo) *rjuo) *xnC3*19) 

CoMMON/BtK?/OMEGA(3, t 9 f 19) , AP(3tl9, 19J ♦BB(3 »iq. 19) ,0»AR(3M9. 19) * 
lD0AR(3 t 19t 19) ,ABAR(3,19»19) «Y(3»l^t 19) »ALPH(3*t9,l9)*BFTA(3*19,19) 

3 /BLK3/XE»PI*H,CU»XK,XM»RCH»PAZ 

a/PLKa/e( 3 , 19) t EZ(3,l9» 1?) »F(3,l9 t 19) ,G(3 t 19) ,70(3*1 9) 
S/RLK5/T*TR»Nao,XNE,W, A (40»90)*DlVlOE 
5/BLK6/LE»Ii*l2,K7»K8 f K9»lM0«IC0 
7/BLK7/EPS»EpE2f AK(tt) ♦Tl»StBETAf SlC t StAP,SiR z 

fl/BUK8/zANtaO) ♦ xNEQ(«o) 

9/BLK13/NU »NF *NT » Ha ( 3) , NS7 ( 3) *MI C j 00 > * MJ C 1 00> ,MK ( lOo > tNALPH 
CoMMON/BLKl5/NFI(100)*NFj(100) f NFK(i00) *NOF 
wplTE (8*10)*K 

10 F0RMATU7H1&NTER AMAT* KK *12) 

DO 115 IsNFfNL 
N$H * nSTCI + D 
N$L * NSTCI-1) 

NS * NsT 1 1) 

NrTART ■ MACl) 

NEND s MA(I*1) « 1 
IF (1*3) 20.15,15 
15 NEND * NSTaRT 
20 < * 0 

DO HO MbNSTARt»NEND 
k = K + 1 

IF U«NF) 35*35*25 
25 Nj * MfcCl * 1) 

NK * M A ( I ) • 1 
J * 0 

DO 30 n*NJ»NK 
J * J ♦ 1 

A (M » N) * XNE*OBAR(I»t *J»K) ♦ **ABAR(I*1*J*K) 

30 Continue 
35 J M 

Do ^5 N*NST ART * NtNO 
J * J + 1 
YY « Y(I,J*K> 

IF(KK # eQ,0) YY = W 
IF (M»N) 40»50,45 

40 A { M * N ) a XNE*OMEGAtI» J»K) ♦ YY*BB(l*J*K) + Ap(l,j* K ) 

GO TO 95 

45 A ( M » N ) a XNE*OMEGA(I,J.K) ♦ 8B(1*J*K)*YY 
Go TO 95 
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44* 

45* 

46* 

47* 

48* 

49* 

$0* 

51* 

52* 

53* 

< 54 * 

55* 

S6* 

S7* 

58* 

g9* 

60 * 

61 * 

62 * 

63* 

64* 

65* 

66 * 

67* 

68 * 

69* 

70* 

71* 

72* 

73* 

74* 

?S* 

76* 

77* 

78* 

79* 

«o* 

M* 

82 * 

83* 

84* 

85* 

86 * 

87* 

88* 

89* 

9ft* 

91* 

92* 

93* 

94* 

95* 

96* 

97* 

98* 

99* 

lftO* 

lot* 


50 SO® * 0*0 

SaL * 0*0 

SOM s o,o 

SA8 s o.O 

SRV * 0*0 

SO® * ft * 0 

SAP = 0»0 

sbe * 0»0 

IF (I*NL) 55»65*85 
55 DO 60 

SO® * SOB * OBAR ( I f J,L) 

SA® * SAB 4 ABaRCI, J t t) 

60 CONTINUE 

65 IF (I»nF) 80»8ft»70 

7ft DO 75 Lt«l,N$L 

SrE s SBE * BETAlItJ,U) 

SOB * sDB ♦ DbarCI,J,U) 

SAL « SAL 4 ALPH(I,J LU 

75 CONTINUE 
80 Do 90 l*1*NS 

IF Cl.EQ.J) GO TU 85 
SOM ■ SOM * QMEGA(I # j,L> 

Yv * YfUJ.l) 

IFIKK.FO.O) YY « W 
SBY e sBy 4 8B(I»J»L)*YY 
85 IF (L.GE.J) GO TO 90 
SaP * SAP 4 AP(X»J,L) 

9o Continue 

A(M*N) **XNE*(SOW*SAl*SOM)»W*SAB'-SbY-XNE** , '*Spf-XNE*(XnE*SDb)-SAP 

95 continue 

IF <I*NL> tOOftlOiUft 
100 N J * Ma(i+1) 

K1K S NJ 4 NST(I+l) * 1 

J ■ 0 

00 105 N*NJ»Nk 
J * J + 1 

A f M» N) e XNE*(AC.Ph(I* 1 * J»K) 4 W*Bt T A ( 1 4 1 » J » K ) 4 xN|r*DB AR (1 4 l * J * K ) 
105 continue 
HO CONTINUE 
115 continue 

write (6» 120 ) divide 

120 FoPMATcteno divide m a trix byipfi2.4) 
if ckk.ne.o) go to 130 

11 9 ma(nF) 

12 * MA(NL) * NST(NL) » l 
IP08MA(K7)4K«»i1 

ICO * mA(K7> *.114 1 
NOO* I2» 1 1 4 1 

WRITE (6, t25)U»n»I?»IRO*lCO»NOQ 

125 FORMAT (1 HO, 3*, 2HLE f 4x,2HIl,4X,2HI2,3X,3HlRQ,3X,3Hlc<MX ,3^00/616 

130 IF (Il.EO.i) GO TO 140 
DO 135 I at 1 » NOO 
DO 135 J* l * NOO 
lL B Il*l*I 
IM* 1 1 w 1 4 J 
A(ItJ)sA(lL»lM) 

135 A(IL»lM)aO, 

140 DO 145 Irl,NOO 
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l()2* 

103* 

lfl4* 

105* 

1{>6* 

107* 

t 08 * 

io9* 

1 io* 
1 1 i* 
U2* 
U3* 
1|«* 
US* 
1 16* 
it 7 * 

ns* 

1 J9* 
1?0* 
1 ? 1 * 
1 ? 2 * 
1 ?3* 
l?a* 
1?5* 
1?8* 
127* 
1 ? 8 * 
1?9* 
1 50* 
151* 
1 52* 
153* 

13 a* 
155* 
136* 
137* 
I38* 
1 39* 
laO* 
Ul* 
142* 
la3* 

laa* 

145* 

U6* 

147* 

148* 

l«9* 

150* 

151* 

152* 

153* 

154* 

155* 

156* 

157* 

158* 

156* 


00 145 J*1»NU0 
145 A ( I , J ) * A(I,J)/DlVIpE 
WRITE (6,150) 

150 PnRMAT(65H0Q(I> * SU M OF Aj.j. E L EmEnTS IN Column 1 OF RATE EquATI 
j ON MATRIX/) 

00 165 1*1, NOO 
OT*0 * 

DO 155 J* 1 »N00 
155 GT»QT*a(J*I) 

WRITE (6, 160)1, OT 
160 F0PMAT(3H Q(l2,3M) *jpEl5,6) 

1 65 CONTINUE 

IF CKK) 170»170»175 
170 WRITE (6,300) 

GO TO 180 
175 WRITE (6,i05) 

180 DO 185 1*1, NOQ 

185 WRITE (6,190) J, (A(I,j) ,J*t, NOO) 

190 FORMAT (1 HO, 13, 2X»lP9El4,tt/(6Xt9El« t a) ) 

WRITE (6,320) 

DO 200 LC*1«N00 
Do 200 l r R*l»NOO 

IF (tC.GT.2, AND,LR,GT,2) 00 TO 200 
CAP(lRftC) “COF aCT (UR,lC»NOO»A) 

WRITE (6,195)LR»U,CaPUR»LC) 

195 F0RMAT(22H COFACTOR oF ELEMENT A* «?I4,2X # 1H*1 PF 15.6) 

200 continue 

pp* o, 

DO 205 I*t,NQO 
205 p P Bpp*CAp(l»I) 

AlAMB*! ,/PP 

A$UM*0, 

CF1*2*/XNE*(SqrT(2 # *pI*(XM/H)*T)*9GRTUK/H) )**3 
JF « 0 

DO 21S J*NF»NL 
JL * N S T(J) 

Do 210 1*1, JL 
JF * JF ♦ 1 

PPG) * (E( J» I) - t(LE,l))*H*CL/(XK*T) 

P2 * 0 ,5*PP0 
EX * EXP (*P2) 

RI(JF) * CF1**( J*LE)*G(J» I)/G(i E,1)*EX 
Rl(JF) * RI(JF)*EX 
ASOM * ASUM ♦ Pi 1 JF) 

210 continue 

215 CONTINUE 

DO 220 1*1* NQQ 

ZAN(I) * AlAMB*CAP(X,I) 

220 CONTINUE 

write (6,225)ALAMb»AsUM 

225 FqRmATOOHOLAmBOAI *ipE1S,6*5X,9HLAmBQAE *,Ei5.6) 

IF (KK) 230»230»235 
230 WRITE (6,310) 

Go TO 240 
235 WRITE (6,315) 

240 Do 250 1*1, NOO 
LX»IWI«1 

XnEQ(I) z Rl(I)/ASUM 
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|U0* 

Ul* 

ih2* 

163* 

16a* 

165* 

166* 

167* 

U8* 
16** 
170* 
171* 
172* 
173* 
17«* 
175* 
176* 
177* 
17«* 
17** 
1B0* 
1 A t * 
182* 
183* 
184* 
185* 
186* 
187* 
188 * 
189* 
1 60* 
161 * 
162 * 
163* 
164* 
165* 
166 * 
167* 
168 * 
169* 


RjU) * zancd/xnegcd 

WrITE ( 6 , 245 ) Lx »ZAN ( T) » R J (I ) *XNEQ ( I ) 

245 F 0 R*AT(l 3 » 2 X» 3 HN * , i pE 1 5 • 7 ♦ 2 * ♦ 7 HN/NEQ * tF- 15 . 7 , 2 X , 5 hNEQ *,Ei 5 * 7 ) 

250 CONTINUE 

GO TO (255,26o),K7 
255 SIG23 a 3 f 74t»l2/SQRT(T) 

Xk|20 s ZAn(1)*XN£/10.0 
GO TO 265 

260 SIG 23 s 2 , 84 t»l 2 / 30 RT(T) 

Xn 20 * ZAN( 20 )*XN£/ 1 o ,0 
265 HH * T l /(XN 20 *Sl& 23 ) 

WRITE ( 6 , 270 ) HH 
270 FoRMAT( 4 H 0 H * 1 PE 11 , 4 ) 

IJ ■ 0 

00 280 I»NF»NL 
JL * NST(I) 

CO ?75 J*1,JL 
IJ * IJ ♦ 1 
XN(IfJ) » ZAN(IJ) 

275 CONTINUE 
280 CONTINUE 

WrITE ( 6 , 285 ) 

285 FoRMATCljHiC I, J, K ) , 5X , IMF, U X , 2H0D) 

DO 295 L»t,NOF 

1 * MF I ( L ) 

J * nfj(u 
K » NFk(I) 

OC * (xN(I»K)/xN(K7,K9)5*(F(I*K,J)/F(K7,K9,K6j)*(E7(K7 f K6,K8)/ 

jE7(I»K,J))*Tl 

write ( 6 , 260 )l,K»JfF(pK,J) ,00 
290 Fo« ma T( 2 H (I 2 ,lH*I 2 HH»l 2 *lH),lP£t 0 . 3 »E 15 . 7 ) 

295 CONTINUE 
return 

300 FoRMAT( 1 HI, 2 X,iHI, 50 x» 20 HRATE equ*tIo n natri*) 

305 F 0 RMAT(lHl f 2 X,iHi, 3 fx*^ 2 HRATE EQUATION MATRI* INCLUDING t*ICK RAD) 
1 AT I ON FIFED) 

310 F 0 RMAT( 24 HlQRTiCALtY THIN SOLUTION/) 

315 FORMAT(44H1SOluTION INCORPORATING THICK RADAilON FiELO/) 

320 FoRNAT(lHl) 
end 


IAoNOSTIcS 
MoN TIME S 


3,87 CPU SECONOS 
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l* 

2 * 

3* 

o* 

5* 

7* 

s* 

9* 

1 0* 

U* 

I 2* 
t** 

10* 

|5* 

1 6 * 
j 7* 

16* 

j9* 

? 0 * 

gnostic® 
ON TIM? 


BLOCK DATA 

C0NM0N/BLK3/Xc,PI»Hi C L*XKi xm»Rch*paz 
fl/BLK4/p(3»l9),eZ<i«l9i I 4 ?) *F (3» 19,19) ,G ( 3 , 19) , 2 R ( 5 , j9) 
9/BLK13/NL*NFfNT»HA(3) ,NST(3) ,MI(I00) .HJ(tOO) , MK ( \ 0q> ,NALPH 
OaTA MA^J»20»30/»NST/I9»i0,i/ 

OAT A CcQ(! v J)"Mtl9),!9tf3)/U«3' 9 l t «9 c9 3M3,«t' 9 Q 1V iS.'5,t3.«3 
1 l»*9,,l5,»S,*27,*7»,3,* 

22, *2»»6»»2**6,, tO* , 2, , 6, , 1 0 • » 1O,*9*0,» 

3t.»18*0*/ 

OATA { ( E( l , J) * ♦ 19) ,X«l,3)/0, # tS9850.32»t 6*271, 7o* t6908t,5o, 

1 171 t29,t5» 18323 1 *08, J84859.06,i65!»59,09 f tb6095t90»i86099*22* 

1 186203,62 1 1 90292.96, t90930, 50 *l9t 2ll* 42, 19^38. 83*1 9io«0.7t, 

U9l 006, 61 *1 91407*20,191 086, 95, 

2 198305, 00 *527460 *57 ,627087, 02 *588045, 76 , 588446.09, ?88407, 6a* 
2609781, 98*609782. 28, 609782, 77, 609702, 95, 9*0,0, 

3637213,67*18*0,0/ 

DATA ((ZB(1*J),J b 1,19)»X b 1,2)/1 9*2, ,10*1., 9*q./ 

Da 7 A XE*PI*H*CI**K*XM,RCH*PaZ/0,803E-10*3.14i6*6,6?56e»27* 

1 2,9979El0.1,3805E«t6.9,lo9tE-28*2,l79E-lt,8,797E-17/ 

End 


,52 CPU SECONDS 


1* 

2* 

3* 

4* 


5* 

6 * 

7* 

8 * 

9* 

i®* 

U* 

12 * 
1 3* 
10 * 
1 5* 

l6* 

i 7 * 

i«* 

1 9* 
20 * 
21* 
? 2 * 
23* 
20 * 
25* 
?6* 
?7* 

?«♦ 


function cofact(Nr,n c ,ne,0) 

DIMENSION 0(Oo,00)*E(00*00) 

MK| * On 

NP»NE«1 

Do to I*1*NP 

IR S I 

IP( I ,GE,NR) IR*IM 
00 10 J* 1 » NP 
IC*J 

IfCJ.CF.NC) IC»J*1 
10 Ef I , J)*D( IR.XC) 

Y ■ 1.0 

M = NDETRM(NN,nP*E,Y5 
IF 30*15,26 

15 WRITE (6,35)NR,NC,Y 
Go TO ?5 

20 WRITE (6,00)NR,NC 


25 Y * 0,0 
30 LT * NR + NC 
MZ«M0D{LT*2) 

Z*t , 

IF(MZ.nE.O) Z»«l* 

Co? act«z*y 

RETURN 

35 FORHAT(23HO****** 
i y »ipeu,o) 

00 FORMAt(20HO****** 
END 


OVER/UNDERFLOW*10X*OHNR *I3*10X,0HNC al3*l 
SINGULAR MATRIX, IOX.OHNr 8i3,lO x *OHNC *13) 
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l* 

2 * 

5* 

a* 

5* 
6 * 
7* 
8 * 
9* 
10 * 
11 * 
1 2 * 


tS* 

tb* 

1 ?* 

t** 

19 * 

20 * 

21 * 

22 * 

?3* 

2«* 

25 * 

26* 

27* 

2 ®* 

? 9 * 

30* 

51* 

52* 

33* 

34* 

35* 

36 * 

37* 

39* 

39* 

40* 

at* 

42* 


FUNCTION COLL (T, J » J»K*N) 

DIMENSION EVC3,l9 f 19) ,A(3f 19»19),AN(3*19, 19) ,*L(3»l9 # l9) f AB(2) i 
1ABN(2),A8L(2> 

CO M MON/BLKJ/XE,Pl f H,CL*XK,XH f RCH,PAZ 

C0^M0N/BLKa/fe(3» 19) * e zf3» 19. 19)*F(3»l9, 19) ,Gf 3»19) # Z8(3*19) 

DATA CACl.l.L),L*2.l9)/U52E*l0*6*68E*H*l,3aE-10*6*0ie*lS. 
l5*03L»l6 f l»O e E*i7 t 3,poE»il « 8 • 1 1 E* 13 « 6 .24E« j 8 , g • 36b- 1 5 , 1 « j 4E*15 . 
?1 .18f :.i7 f 5.2bE # t2,2,25E«l3* U®7£«2ti3«0lE*l5 t l ,36E. 18,4,88 f.* 16/ 
DATA (A(l*2*L) f L*3fl9)/l • t0E<*0fe»9.47E<*Q9,4,5(,E»O8,4.3aE«.09, 
l6.O3E^o9.S,3«E-O9,0.36E*o9 f 7.72E«O7»4fO6b«*o7 t 7.27E-iO f l*O®E>O9t 
?3«99fc<» 1 0 . 1 ,2lE*09,8«44E»O8,2*46E«’l0, 1 •^®E»08, 3. 75E«*08/ 

DATA (AC1.3.L) ,L»4,5)/4.45E*07 f 1.00E-07/ 

DATA C AN( 1 » 1 *L) » L«2» ) 9) /O • 927*0 *9 9 8?0««76tl, 39* i,4i» 1.68. 0,989. 
t'O. 520. 1,72, 1 . lb. 1,29 .1,61 *0.594, 0*580 *2. 32, 0,555, 1,27, 1.39/ 

DATA (AN(1 ,2.L) ,L»3»i 9) /»0. 196, 0.586, 0.154, 0,390. 0.348, 0,287. 
i0,394 V( ,0,164.»0, 135*0,41 3, 0.230.0.410* 0.441, „o,o5l4.0,423* 

2« 0*0663**0. 024/ 

DATA CAN(1,3*L) »L*4 *$)/•• 0617, 0.389/ 

DATA (AL(l*l*L)*L»2tl9)/0.999,l, 0*1.0 1*1.0,1.01 .0.995,1,01,1.02, 
1 0 , 997 .1.0*1*01, 0.997, 1.01,1.02, 0.976.1.02,1, oi* 1,0/ 

DATA ( AL (l*2*L),L«3»l9)/1.4l, 1,36* t,28. 1.15, 3*1, 12,2*1, 18. 1,07, 

1 1.14*1.05*1.12*1 *09,) ,13. 1,09*1.10/ 

DATA (AL(1»3*L)*L*4,5)/1.43*1.72/ 

DATA caB(L) ,L*i,2)/6,99E*1S*1,02EM4/ 

DATA (ABN(L)*L«i*2)/i, 44,1,23/ 

DATA (aBl(l)*l«1*2)/1, 0*1,01/ 

XKT * XK*T/1.602 £ "12 

EV(I*J,K) « i,23977E^04*EZ(I,J,K) 

Go TO (10*15),N 

10 EX * At (I.J,K)*EV(I, j,K)/XKT 
At * T**AN(I*J,K) 

A2 a ExP (*EX) 

COLL * A(I,J,K)*Ai*A? 
return 

15 EF * A8S(E(I+1,K) • £ ( I *J) ) 

EG = 1,23977E* 0 4*EF 
EX * A8LtI)*EG/XKT 
A) a t**ABN(I) 

A? a ExP(»EX) 

CqLL « A0(I)*A1*A2 
RETURN 

end 


IAgNOSTIcS 
TION TINE * 


,90 CPU SECONDS 
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1 * SUBROUTINE C0N3TS 

2 * COMMON /BCK3/XE*PI*H,CI»XK*XM,RCH*P4Z 

34 3/BlK2/0MEGA(3,t9U9) , AP (3* 19* 19) »BB (3* 19 , 19) ,OBAR ( 3 , 19 ♦ 19) ♦ 

4* aDB*Rt3»19tl9) f A 0AR(3, 19,19) *Y (3*1 V,19J»ALPM(3, 19*1 9 ) »BET A (3*19,19) 

5* a/BlK4/E(3*l9),EZCl*i9,19)*F(3,l9,i9),G(3,19), Z BC3*l9) 

69 5/BLKS/T*TR*NOO t XNE*N,A<<10*40),OlVX0E 

7* 9/BLKl3 /NL*NF ,Nt*MA(3) ,NST(3)*MIC100) ,MJ ( 1 00 ) ,MK ( 1 0$) *NAIPH 

8 * DATA B, AX*0*OP,EP* ADp*DDP/^, 0625* *932» *Oh25*.o04* 17**.15**025/ 

9* 10 foRMAT(1hO,12x,3m1 «,2X* 1H2* 19X* 1M3//«X* 1HK*8X*3hJ »* 2 X * 1 H 1 * 1 9X * 1H 

io* it/) 

119 15 FORMAT(I5»1P2E20,7) 

129 20 FoRMAT(1hO, 12X,3HI «, 2X . 1 HI * 19X * 1H2//4X * 1 H J* B x ♦ 3HK «,2X, 1 Hi * 19X . IP 

139 ll/J 

*49 25 FORMAT( 6 HOX * J) 

1 54 30 FORMAT(lH0*3X,lHa*6X,3HK ■ * 2X , 1 HI * l 7X * 1 H2* 1 7* , 1 H3* 1 7X , 1 H4 * l7X * 1H5* 

164 1 17X*1H6, 17X* 1H7/3 

j74 35 FoRHAT(1H0,10X*3HK « , 2 X , 1 H 8 . 17X * 1 H9 , 1 7X , 2H 1 0 , 1 6 X *2h1 1 * 1 6 X *2H12, 16X 

1 84 i*2Hl3,l6X*2HU/) 

j94 40 FORMAT(1HO*3X*1HJ*6X,3HK ■ * 2X * 2Hl B* 16X , 2H16 » 1 fcX * 2 H) 7* \b* ,2Hl8, 16 X • 

209 1 

214 45 

22* SO 

239 1 

244 55 

254 60 

264 65 

274 70 

284 
294 
304 

31* 

32* 

33* 

34* 

35* 

36* 

37* c 
38* C 
39* 

409 

at* 

42* 75 

43* 


2H19/) 

F0RMAT(6H0I * 2) 

FORMAT (IH0*3X»iHJ, 8X,3HK « *2X * t Hi * i9X ♦ 1 H 2 * l9* , 1«3 * j 9X , 1 M 4 » j 9X , iH5 * 
/) 

FORMAT (1H0,12X,3HK « , 2X , 1 H6 * 19X * 1 H7 * 1 9X , 1 H8 , 1 9 X , l Hg* l9X,2H t 0/) 
FoRMAT(I5*lP7El8,7) 

FORMAT(X5*lP5Et8.7> 

FORMAT(I5*1P5E20,7) 

TE » T 

TpTE b TR/TE 

FF» * CSQRTC2,4PI*CXm/M)4TE)*SQRT(XK/H))**3 

FF2 * FFl*TRTE**t.5 

HCKT*H*CI./(X«4TE) 

HCKTR e HCKT/TRTt 
TA» = ABStTRtE - 1,0) 

Ti*8QRT(TE) 

fff*i,/ffi 

4*44***49*4**449 ORlGlNAtLY SUOBA 9*444449********4**44***9*4**< 

DO 130 IJ*!*NAUPM 
IF (TAB.LT,1.0E*05) GO TO 80 
READ (5»75)I*K,J.TTR f AL*BE 
F0RMATC3l6 f 3fcl5,8) 

BETA(I*K,J) * BE 
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46* 

a7* 


48 * 

49 * 

50* 

51* 

52* 

53* 

5«* 

55* 

56 * 

57* 

Sft* 

S9* 

60 * 

61 * 

62 * 

63* 

64* 

65* 

66 * 

67* 

68 * 

69* 


70* 

71* 

72* 

73* 

74* 


75* 
76* 
77* 
78* 
79* 
80* 
8 1* 
82 * 
63* 


84* 

85* 


86 * 

87* 

88 * 

89* 

90* 

91* 

92* 

93* 


94* 

95* 

96* 

97* 

08* 

09 * 

loo* 

loi* 


GO TO 85 
80 I * MI(IJ) 

K * MJ(IJ) 

J * HK f IJ 3 
TTR ■ TR 

A|_ * ALPH(I>K,J) 

Bp a B£T A ( I , K , j) 

85 I « I - 1 

EF»ABS(E(J*i»K)«Ml*J)) 

X a J 

E?2»EF*CL*H/l.602iE»l2 

X2r*EP*HCKT 

FF3 « At ♦ Bfc 

ARG a pF*HCKTR 

IF (ARG*GE.ae,o) Go TO *0 

Y2 a exP(»ARG) 

GO TO 95 
90 Y2 » 0,0 

95 A0AR ( I , J » K ) * 2,0*G(I*1,K)/G(I,J)*Y2*FF2*FF3 
IF (XZ,GE 9 88,) GO TO 100 
Yt*EXp C-XZ) 

GO TO JOS 
100 Y1«0, 

105 F 1 a (5,1 » 1,2/X «, # 9/X**2) 

IF (CJ,Eft.i)«AKiD*(K t FQ,n) GO TO ItO 
OBAR(I f J f K)al,l5U8*pi*Z8(I»J)*Tl*Yl/tz2**2 
GO TO 115 

110 08AR(i,j,K) * C0LKT,i»J,K,2) 

115 CONTINUE 

IF (XZ,GE,88.) GO TO 120 

DBAR(I*l»K.J)a 6 (I f j)/ G (I+lfK)/2.*bXP(XZ)*FFF*0BARll*J»K) 

GO TO 125 

120 DBAR(I*l*KtJ)*o,0 
125 continue 
130 CONTINUE 

write (6tl35) 

135 F0RMAT(12H1BETA(I» J»K)*27X,42H8TIMUCATE0 RECOMBINATION COEFFICIENT 

i at trj 
write (6»uo) 

140 FORMAT(1HO, 12X.3HI a, 2X, 1H2« 1<>X» 1H3//4X , IHK.gx ,3HJ a, 2 x, jHj , 19X, lH 
1 l/) 

WrItE (6« U5)(K*BETA(2f ttK)'BETA<l'itK)'KM 9l 9) 

145 FORHAT(I5»1P2E20«7) 

WRITE (6 1 150 3 

150 F0RHAT(l2H10HAR(ltJ* K ) »27X # 43HElE C TR0N COlUStONAl IONIZaTjON COEF 

ificient) 

WrITF (6,20) 

WRITE (6,15) lJtO«AR<i,J f l) f oBAR(2#J,l),Jal,l9) 

WRITE (6,155) 

155 FORMAT(12HlOBARa* J, k) ♦ 27X*46HELEGTR0N COtLlSjONAl RECOMBINATION C 
COEFFICIENT) 

WRITE (6,10) 

WRITE (6,15) CK,DBAR(?,l,K),0BAR(3*l,K)fKal,l9) 


WRITE (6,160)TTR 

160 FnRMAT(12HlABAR(I*J,K) #24X,33HRADIATIVE IONI^aTION RATE AT TR slPE 
111.4) 

WRITE (6,20) 


C 
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10 2* 
103* 
104* 
1 Cl5* 

I ft6* 

io 7 * 

100 * 

iO<»* 

110* 

111* 

iia* 

H3* 

H4* 

ll5* 

116 * 

II 7* 
MS* 
1|9* 
1<?0* 
1?1* 
12a* 

J?3* 
i?«* 
125* 
1?6* 
1?7* 
128* 
1?9* 
1 30* 
131* 
132* 
133* 
134* 
135* 
13** 
13 7 * 
138* 
13** 
1 40* 
141* 

lq2* 

143* 
1 44* 
145* 
148* 
147* 
148* 

1 49* 
ISO* 
t S 1 * 
132* 
153 * 
104* 
155* 
156 * 
107 * 
13 «* 
13 ** 


C 

C 


C 

C 


c 

c 


WRITE (6*15) Cj*ABAR(i ,J»l)*A0AR(2»J,l)*J*l,l9) 

M************** ORIGINALLY COMB0 I***************************** 
Xk2 s XK*SGRT(XK) 

T2 ■ T*S 0 RT(T) 



C38a32,*Pl/8QRT(3,)*pAZ/8QRT(2,*Pl*XM)*RCH**p/XK2 

C39a4 t *RCH*PAZ/(S0RT(XM*TE)*SQRT(2 t *pi*XK)) 

WrItE (6,165) 


165 F0«MAT ( lHl*3X # 1 Hl,4X f tHJ,4X,lHK, liX, 8 HB(I,J f K) f l 2 X,0HB(I,K # J)/) 
Dp 200 1*1*3 
JL * NST(I) 

DO 280 J> 1 ,JL 
00 280 Kb 1 « JL 

Ev * l # 23977 Ew 04 *EZ(TfJ*K) 

ENZ * EV/13,6 


XZ*EZ(I*J*K)*HCKT 
Ext * EXP(-XZ) 

EX2 * FX1**2 

Ex 1 0 * EXl**l0 

Ex 17 « EX1**17 

EX30 = EX1**30 

EXIOO b EX10**t0 

IF (XZ.GE# 88 « 0 ) GO To 170 

Y i*EXP («XZ) 

GO TO 1 75 

170 Y 1 *0 • 0 

175 IF (J.GE.K) GO TO 200 

IF (I.GT.l) GO TO 215 

IF (J.EQ.l) Go TO 245 

IF CJ.EQ.2) GO TO 245 

IF C(J,EG.3)*AND*CK*EQ*4)) GO TO 245 
IF (CJ,Eg.3)*AND*(K*fg*5)> GO tO 245 
IP (F(I*K*J)*LE.U0Ep7> GO TO 1 9 0 


*********** ORIGINALLY 8U36 ************************************** 
Y2 s Y 1 /xZ 

z«t./xz 

IF (EP*XZ.GE.08,O) Go TO 100 
Y 3 *EXP(*XZ*EP)/XZ 
go TO 105 
180 Y 3*0 1 0 
185 Fl * 0,t4*EXt 
F2 8 **0*083*EX2 
F3 * »0*053*Exi7 

OMEGAClf C38*F(I» J,K)*(Fl ♦ F 2 ♦ F3) / ( t2*XZ**2) 

GO TO 250 


********************* ORIGINALLY 
1 90 CbC39/gU*J) 

Yo«RCH/(XK*TE) 

XZl*XZ+YO 

IF (XZ1.6E.80O 60 To 1*5 
F1«EXP(»XZ1) 

F j*F3* ( l ,*Y0) 

GO TO 200 
195 F 1*0*0 
Fj*0 • 0 

200 tw CXZ.GE.08O GO TO 205 


SU37 *************************** 


i 

{ 

I 

I 
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> 0 * 



F 2 *EXP(»XZ) 

( 

*f 6 l* 



Go TO 2 l 0 

( 

162 * 


205 

F 2 a 0 • 0 

( 

163* 


210 

OMEGA (i,J,K)»C*((F2*Fl)/Y0*F3) 

( 

164* 



OmEGA(I*J»K) * OHEGA(I» J»K)/ENZ 

( 

165* 



GO TO 250 

C 

166 * 


215 

IF (F(i, K, J) »LE, 1 ,E«t) GO TO 240 

( 

U7* 

c 



t 

168 * 

C 


^* ********* OrX^InAlLY SUS 8 ************************* ************< 

169* 



Y 2 *n/xZ 

( 

170 * 



IF ( 2 **XZ»GE» 00 tO) Go TO 220 

( 

171* 



Y3*EXP(»2**XZ)/XZ 

( 

172* 



GO TO 225 

( 

173* 


220 

Y3°0* 0 

( 

174* 


225 

IF C48.*xZ,GE,88.) Go TO 230 

( 

175* 



Ytt»EXP(«.a 8 ,*XZ)/XZ 

( 

176* 



Go TO 235 

( 

177* 


230 

Y4*0* 0 

( 

178* 


235 

Ft a 0.2*EX1 

( 

179* 



F? ■ (0.045/XZ)*EX2 

( 

ISO* 



F 3 * .(0,024/XZ)*£X1o 

l 

181 * 



F« * *(0.016/XZ)*EX3o 

( 

182 * 



F5 « *(1.33 ♦ o*0O5/xZ)*EX1OO 

( 

183* 



OmEGAU»J»K) a C48 *F(I*J»K)*(Fi ♦ F2 ♦ F3 * F4 ♦ F5)/(t2*XZ) 

c 

184* 



GO TO 250 

c 

185* 

c 



( 

1 R 6 * 

c 


*********** ORIGInHEY SU39 ***************♦»************♦*******( 

187* 


240 

0HEGA(l,J,K)»C39*Y l / e Cl l J) 

c 

188* 



OMEGA(I*J*K) * OMEGArl*J*K)'ENZ 

( 

189* 



GO TO 250 

( 

190* 


245 

OmEGA(i»J»K) a C0 U.(t»I* J*K» 1) 

c 

191* 


250 

IF (XZ.GE#88.) GO TO 255 

c 

192* 



OmEGAu,K*J)*G(I»J)/G(I»K)*ExP(XZ)*OM£GA(I,J,k) 

( 

193* 



GO TO 260 

( 

194* 


255 

OmEGA ( I»K tJ)*O t O 

( 

195* 


260 

IF (F(1,K»J).LE.UE«7J GO t° 280 

( 

196* 



xmu»cl*ez ( i*k* j) 

c 

197* 



CaPb8,bp1**2*{xE**2/xM)*(XNU/CL**2)*(XNU/CU 

c 

198* 



AP(I,K,J)«GU, J)/G(I X)*CAP*F(l f Kf J) 

c 

199* 



BVUK a 4*0*Pl**2*XE**2*Ftl»K* J)/l CXM*CU*CH*XNU)) 

c 

200* 



8 V 1KJ a GCI#J)/G(I,K)*BvIJK 

c 

201* 



WRITE (6,265) I,J»K,BvlJK f 0VIKJ 

( 

202* 


265 

FoRmAT(3i5,1P2E20.«) 

< 

203* 



ArG a H*XNU/(XK*TR) 

< 

204* 



ArG a o,5*ARG 


205* 



IF (ARG-68,0) 275,270,270 


206* 


270 

88(IfJ,K) a O.o 


207* 



B8(I»K,J) a 0,0 


2 ()8* 



GO TO 280 


209* 


275 

PART a EXP(ARG) 


2l'0* 



P L FT a (2tO*H*XNU)*(XNU/CU)**2/(pART * 1.0) 


l 2 i** 



8B(!*J,K> * (PLPT*BV|JK)/(PART ♦ 1,0) 


2l2* 



B0<XtK,J) a (PLFt*8VlKJ)/CPART ♦ 1,0) 


2 j 3* 


280 

continue 


2i4* 



WRITE (6,285) 


2l5* 


285 

FoRmAT( 13H IOMEGA II, j t K) »26Xt 43HE|.ECTR0N COLLISION^ ExCItAtION 

CqE 

2j 6* 



fFFlCIENT) 


2 1 7* 



WRITE (6,25) 
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2l8* WRITE (6*30) 

219* WRITE (6f60) (J»(0MEGa(1* J«K>«Ks1'7)'JBW|9) 

2?0* WRITE (6,35) 

221* WRITE <6, 60) CJ,<UMEGA(l*J*K),Ka8,l4), 4*1,19) 

2?2* wrITE (6,285) 

??3* WrITE (6,25) 

2g4* WRITE (6,4 q) 

225* WRITE (6,65)(J,(0MEGA(l»J*K)fK*i5M < ))»J=l»l9) 

226* WrITE (6,45) 

227* WRITE (6,50) 

228* WRITE (6,70)(J,(OMEGA(2tJ»K)tK«l,5),J*l,10) 

229* WrITE (6,55) 

230* WRITE (6,7o) CJ, (UMEG a(2iJ»K) ,K« 6, iO),J*l*10) 

23 i * WrITE (6,290) 

232* 290 FO«WAT(!0HIAP(l»J»K),29Xf43HEIWSTElM SPONTANEOUS TRANSITION PROBAB 

233* tUlTY) 

234* WRITE (6,25) 

23?* WRITE (6,30) 

236*- WRITE (6»60)(J,(Ap(l,J,K),K»i,7),Jsl,19) 

237* WRITE (6,35) 

238* WRITE (6,60) ('JfUPIl, J,K) ,k» 0»1®) *J*itl9) 

239* wrITE (6,290) 

2« 0* WrITE (6,25) 

241* WRITE (6,40) 

2fl2* WRITE (6,65)(J,(AP(1 ,j,IO*K*15,19),J*1,19) 

243* WrITE (6,45) 

244* WRITE (6,50) 

245* WRITE (6t70> l J»(AP(2,J*K)tKM»5)»4 !B l,t0) 

246 * WRITE (6,55) 

247* WRITE (6,70)(J,(AP(2,J,K) f K*6,i0)»J»l*10) 

248* WRITE (6,295) 

2«9* 295 F O RHAT(10NtBB(i,J,K) # 29X,42HElNSTElN ABSoRPTjqN TRANSIT^ PRoBABI 

250* iLlTY) 

251* WRITE (6,25) 

252* WRITE (6,30) 

253* WRITE (6,60) (J,(BB(1, J,K),K*l*7)t4*i»19) 

254 * write (6,35) 

255* WRITE (6 t 60)U,(bB(l # J,K),K«8,l4)tJ«lil9) 

256* WRITE (6*295) 

257* WrITE (6,25) 

258* WRITE (6,40) 

259* wrITE (6,65HJ,(B8n t J,K),K*15 # 19),J«l,19) 

260* WRITE (6,45) 

261* WrITE (6,50) 

262* WRITE (6,70) (4,(B8(2,J,K),K«1,5),4»i,10) 

263 * WRITE (6,55) 

264* WRITE (6,70)(J,(B8(2 t J,K),KR6,i0)f JaitlO) 

265* Rf turn 

266* end 


UAcNOSTIrS 
V T I0N TIME 
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1* SUBROUTINE EUlM(HHOtAP) 

?* c *****#tHI$ SUBROUTINE COMPUTES DOUBLE tOFACTORS AND USEs THEM 

i * c ***** *to compute eta and jota 

4* COMMON/BLK5/T f TR»NOO,XNE f H f A(ao»40) ,OlVJOE 

5* 1/BUK6/lE»U »I2*K7,K8,KR»IROf ICO 

6* ? /BLK7/ EPS»El ,E2»XK( 4) , Tl *BE tA ,C , Aq*KZ 

7* DIMENSION RMAT(40,40)*IO(4) 

0* DATA lD/lHu»iHvt2HS4,?HS5/ 

9* Mn » 40 

jO* N002 = NOO • 2 

1 1 * SaO, 

12* Sl*0, 

13* WRITE (6,10) 

i 4* io Format (S4Hodoubl e cofactors used To compute eta and Iota In ElIm/6 

,5* tlH SIGN FOR DEFINITION OF COFACTO# IS (•U*»(Rl ♦ p2 ♦ Cl ♦ C2)/48 

t6* ? h this sign changes whenever r 2 > ri or q ? > ci//ion eliminate*? 

1 7* TX.UhRI Cl r 2 C2*21 X , 1 1HR 1 Ci R2 C2> 
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i«* 

19 * 

? 0 * 

? 1 * 

? 2 * 

?3* 

?4* 

?5* 

26* 

? 7 * 

? 8 * 

? 9» 

3®* 

31* 

32* 

33* 

SO* 

35* 

36 * 

37* 

30* 

3 9« 

00* 

ol* 

« 2 * 

o3* 

a«* 

o5* 

06* 

«7* 

08* 

a9* 

so* 

51* 

52* 

S3* 

59* 

?5# 

50* 

57* 

50* 

59* 
60 * 
bi * 
62* 
63* 
64* 
65* 
66 * 
67* 
60 * 
60 * 
70* 
71* 
72* 
73* 
74* 
75* 


DO 60 KM«1»N00 

IF ((KM,EQ,1R0)*0R,(KH*EQ.IC0)) go to 60 
C******gUBRO!|TlNE ELM2 ELIMINATES 2 ROwS AND 2 COLuMNg 
C******fROM ORIGINAL MATRIX AND PLACES Th£ N*2 MaTri* jNtO 

c******pmat and tsen the Double cofactor is computed 

CALL ElM2(lR0 f tC0 f KM,I«0,RMAT,Z) 

U * 1,0 

u * NDETRM(MN,nOU2»RmAT*U) 

IF («*»2) 30 *15,20 
15 WRITE (6»H5)Ir0*IC0*IDC1) 

Go TO 25 

20 WRITE (6,t20)lR0,JC0,lDCt) 

25 U a 0,0 
30 U * U*Z 

CALL ELM 2 (ICa t lRO # KM,lCO f RMAT t Z) 

V * 1,0 

M » NDETRM(MN f N002,R M AT»V) 

IF <M»2) 50*35,40 
35 WRITE (6,U5)!R0*IC0 t lD(2) 

GO TO 45 

40 WRITE (6, 120) I RQ* ICO, ID (2) 

45 V * 0.0 
50 V » V*z 

IF(KM.GT.IRO) U * »U 
IF(IR0,GT,IC0) U » «u 
1F(KM,GT,IC0) V * *v 
IF(ICO,GT,IRO) V a 

C ****** A IS THE MATRIX ELEMENT, U ANp V ARE THE DQubLE CoFAcTOrS 
c *M***tHE MATRIX WAS divided by DIVIDE to PREVENT OVERFLOW when 
c ******the cofactors aRe being computed, this must be t 4 ken into 
c ****,*aCCOUnT here 

S*S+A(kM,IRO)*U*DIVIdE 

SI»S1+a(KM,IC0)*V*DIv1DE 

WRITE (6,55) IRO* ICO,kM»IRO,UiICO»IRO*KM, IC o*v 
55 Format ( iox, 413, 2X, IPE15. 6,5X* 413, 2XfE15, 6) 

60 continue 

call ELH2(IR0,ICU,ICo*IR0*RMAT,Z) 

sa * 1*0 

M » NDETRM(MN,nOQ2,RmAT,$4) 

If (M*?) 60*65,70 
65 WRITE (6,115)Irq»IC0,ID(3) 

Go TO 75 

70 WRITE (6,U0)lR0*IC0 t !D<3) 

75 SO » 0,0 
60 SO » S4*Z 
SO « -S4 

write (6, 105) IRO* ICO, ICO* IRO, So 
Call ELM2(ICO,IRO,IR 0 ,ICO*RMAT,Z) 

$5 * 1,0 

M B NDETRMCMN,n002,RmAT,S5) 

IF <M»?) 100*65*90 
65 write C6,US)IP0*ICQ,IDC4) 

GO TO 95 

90 WRITE C6,120)XRG*ICO,ID(4) 

95 S5 ■ 0,0 
100 S5 a S5*Z 

SS * m §5 

WRITE (6,10S)ICO*IRO,IRO,ICO*S5 
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^ 6 * 105 FdRMAT(10Xt«l3f2X*lPEl5.6) 

>Jf 7* ET**1 */Sfl*S/Ap 

78 * XlOTA»RHO/S*l*Sl/AP 

79* wrItE ( 6 « 1 10) Et A * X IOf A » EPS 

80* 110 F0RMAT(iH0#5HETA *lPn5*6* | OX »6H20TA *E 1 5*6* t ©X ,sHfrPS sEj5,6) 

81* WBlTECMlU) S»St 

82* 111 FoR ma T(<*H S »lpEU,7 t lOX.«H81 *El«,7) 

83* WRITE C6,U)Tl,TfW|Ei,E2»TR 

8«* 1 i F0RMAT(7x»2Htl t l«X»2HTE*14X#lHWf t^Xt2HEl*U)( t 2HE2»i«X,2HTR/ClP6El6 

85* t. 6 )) 

Ah* RETURN 

«7* 115 FORMAT(2lHO******UNOER/OVERFLOWtl«X»5HROW *13* tOX* 5 HCOl «I3*10XtA6 

88 * 1 ) 

89* 120 FoRMAT ( a 2 Ho**»*** 8 IN(;ULAR MATRjX, loX ,5HRUW = 13 , 1 0 X,5HCOL «!3 t |0X*A 

90* 16 ) 

91* END 

^gnostics 

Ion time • 1.R6 CPU SECONDS 
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1* 


SUBROUTINE ELM2(lRA»TCA*KRA*KCA»RMAT*Z) 

2* 


CoMHON/BLK5/T,TR*NOO,XNE # W f A(40 i«0) , DIVIDE 

3* 


DIMENSION RMAT (40*40) 

a* 


IR0*IRA 

3* 


IC0»ICA 

6* 


kpObkra 

7* 


KCO*KCA 

e* 


IF (IR0.IE.KRQ) Go To 10 

9* 


I$R*IRO 

lO* 


ISC«KR0 

1 1* 


IrOsISC 

12* 


kro«isr 

! 3* 

10 

IF (IC0.LE.KC0) 00 To 15 

14* 


isr*ico 

15* 


iaC«KCo 

!*•* 


ico«isc 

17* 


KcO»ISR 

{8* 

15 

NN*N00-2 

j9* 


DO 45 IJK«I»Nn 

20* 


IF (IJK.GE.IRO) GO To 20 

?1* 


NP»I JK 

32* 


GO TO SO 

?3* 

20 

IF (IJK.GE.KRO-IJ GO TO 25 

24* 


NP»JJK+1 

?5* 


GO TO 30 

2 6* 

25 

NP«IJK+2 

?7* 

30 

Do 45 KO«l»NN 

28 * 


IP (KO.GE.ICO) GO TO 35 

29* 


NR*KO 

JO* 


Go TO 45 

31* 

35 

IF (KO.GE.KCO.i) GO yO 40 

32* 


NRPKO+i 
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0 3* 60 TO as 

a* 40 NR*K0*2 

35 * 45 RMATCIjK»KO) a A(NP*NR) 

3 6* LT elRA+lCA + KRA+KCA 

37* LZ*H 00 (LT» 2 ) 

38* Z*l* 

39* IfCUZ.nE.O) Z*«l* 

<40* RETURN 

41 * END 

iagnostjcS 

T ION TIME * **>7 CPU SECONDS 
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1* 

function ndetrmcn* ln, a *dj 

2* 

H»0 

5* 

CALL stmeq(n*ln*o,a*o*d*o*m) 

4* 

NDETRMbM 

5* 

return 

6* 

end 

iAfi^OSTlcS 

TJON TIME « 

•21 CPU SECONDS 


1* 

function nsjmeqin^n^lm* a* b*d* e) 

2* 

M*1 

3* 

CaU SlMEQ(N*LN*lM»A*BtD*E t M) 

a* 

MSIhEq»m 

5* 

RETURN 

6* 

end 

AGNOSTICS 

ION TIME » 

,21 CPU SECONUS 


1* 


SUBROUTINE SlMEQ<NitN*LM»A,B*D*E.M> 

2* 


Integer e 

3* 


EQUIVALENCE ($AVt « ISaVE) 

4* 


DIMENSION E(t,N) * A(n*n) »B( NtLM) 

s* 


IF (M.EQ.O) G 0 TO 15 

6# 


DO 10 I«1*LN 

7* 

10 

E(I)«I 

8* 

15 

LnH1«Ln«1 

9* 


call OvERFL £IBIC) 

1 0* 


DO 80 k*1*LNN1 

11* 


SAVE«*i ,0 

12* 


Kt=KM 

13* 


DO 25 J»K»LN 

j 4* 


Do 25 I«K»LN 

15* 


IF (SAvEi,ABS(A(I»J))) 20f 25*25 

16* 

20 

5aVE*ABS<A(I*J)) 

17* 


JBIG* I 

J8* 


JBIGb J 

19* 

25 

continue 

20* 


IF (K*EQ.IBlG) go to 40 

21* 


D*»D 

22* 


Do 30 J*K»LN 

23* 


SAVE«A(K,J) 

24* 


A(K,J)bA(IBIG,j) 

25* 

30 

A<lBlG f J)*SAVE 

26* 


IF (M.EQ.O) Go TO 40 

27* 


DO 35 J*1»LH 

? 8* 


SaVE*B(K,J) 

? 9* 


6(Kf J)«B( IBIG,J) 
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32* 

33* 

34* 

3?* 

36* 

37* 

38* 

39* 

40* 

41* 

42 * 

43* 

44* 

45 * 

46* 

47* 

48 * 

49* 


52* 

53* 

54* 

55* 

56* 

57* 

5®* 

59* 

AO* 

61* 

62* 

63* 

64* 

65* 

66 * 

67* 

6 ®* 

69* 

70* 

71* 

72* 


73* 


74* 

75* 

76* 

77* 

78* 




74* 

80 * 

1* 

82 * 

83 * 

84 * 


35 8 ( 1BIG, J) *SAVE 
40 IF (KwjBjG) 45,55,45 
45 D»«D 

00 50 I ■ t ♦ L N 
$AVE*A(I,K) 

A(I,K)«A(I«J6lG) 

50 A ( I , JBIG) *5 A VE 

IF CM,EQ,0) GO TO 55 
ISAVE*E(K> 

E(K)«E(JBIG) 

E(JBIG)*ISAVE 

55 IF (ACK*K>) 6o,130»6o 
60 0»D*A(K*K> 

00 80 !*Kl,LN 
SaVE«A(I.K)/A(k*K> 

00 65 JsKtfLN 

65 A(I»J)sA(I,JJ»SAV£*A(K»J) 

CALL OVERFL(IBIG) 

IF (IBlG-n 7o,l25,7o 
70 IF (H.EQ.O) GO TO 80 
D 0 75 J«l»LM 

75 BtI,J)«atI,J>»3AVE*B(K,J) 
CALL overflcIbig) 

IF (IBIG.I) 80,l25,8o 

00 continue 

IF CAClN,LN)) 85*130,85 
85 IF(LN*nE# 1) O*0*A(i r N t LN) 
CALL OvERFL(IBIG) 

IF (IBlG-l) 9o,l25,9 0 
90 IF CM) 95 i 120, 95 
95 Do 110 J*l.LM 

BCLN, J)*B(LN,J)/A(LN f LN) 
call overflcibig) 

IF (IBlC-n 100*125,100 
100 00 110 JBIGM,LNMi 
I»LN*JbIG 
3AVE«0, 

1 p l a I * 1 

00 105 K*IP1*LN 
105 8AVM$AVE*A<I»K)*BCK,J> 

B(X*4)«CB(I*J)*3AVI)/AC1*I) 
CALL OVEPFLCXBIG) 

IF (IBlG*l) 110,125,110 
110 CONTINUE 

Do 115 K«1,LN 
I«E(K) 

DO H5 JaltLM 

115 A(I , J)«B (K,J) 

120 M«1 
return 
125 M»2 
return 
130 M«3 

RETURN 

END 


agnostics 
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Input Quantities 


NF 

NT 

NALPH 

NOF 

F(I,J,K) 

T 

TR 

ALPH 

BETA 

Tl 

W 

XNE 

K 7 ) 

k8 ) 

K9 ) 

LAST 

LE 

NTAU 

TAU 

DIVIDE 

NOY 

IPROB 

HAST 


initial ion (usually l) 
total number of ions (usually 3) 
number of a's and p’s 
number of f numbers 
f numbers 

electron temperature (°K) 
radiation temperature (° K) 

a ijk 

^i<jk 

optical thickness at center of line (K7, K8, K9) 
dilution factor 
number of electrons 

indeces of thick line 

not used 

some ion between NF and NT (usually 2) 
number of depth points, TAU 

depth points at which optically thick solutions are made 
constant divisor to prevent overflow of A matrix 
number of thick lines other than line (K7, k 8, K9) 
problem identification number 

= 1 , go to beginning of program and start new problem with new 
temperatures, etc. 

= 2 , start in middle; temperatures remain the same, but new 
T-, W, XNE, etc. 

= 3 > change the Y values only. 


11 
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Input Cards 


Card Format 


Content 


1 1216 

2, 2a, 2b, etc. 316, E12.8 

3 6E12 . 8 

k, 4a, etc. 316, 15X, 2E15.8 


5 

6 

7, 7a, etc. 

8 

9 

10 
11 


6E12.8 

1216 

6 ei2 .S 

6ei2.8 

1216 

1216 

316, E12 . 8 


NF, NT, NA.LPH, N0F 

K, J, K, F. a set of cards containing 

3 3 ijk 

one f number and its indeces per card. 
T, TR 

I. J, K, a. , 0. „ a set of cards 

3 3 3 ljk ijk 

containing d and P for a given ijk for 
a given temperature. 


T x , W, XNE 

K7, K8, K9, LAST, LE, NTAU 
TAU 

DIVIDE 

N0Y 

IERjfe, HAST 

(if needed) Y. ., = 1.0 a set of cards 

1JK 

equal in number to N0Y containing the 
indeces for thick lines and Y„^ = 
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Some Other Constants and Variables 


TE 

AL, ALPHA 
BE, BET 
MA 


NET 

NX 

G(I,J) 

E(X,J) 

ZB(I,J) 

EX(l,J,K) 

ev(i,j,k) 

Bp 

BTR 

BpBTR 

XE 

PI 

H 

CL 

XK 

KM 

RCH 

PAZ 

A(I,J) 

NFI j 
NFJ ) 

NFK ^ 


sometimes used for T (electron temperature °K) 
sometimes used for ALFH 
sometimes used for BETA 

array containing the serial number of the first level in each ion. 

That is, MA l = 1, MAg = 20, MA^ = 30. 

array containing the number of levels in each ion 

number of the last ion (=3) 

statistical weight ion i, level j 

energy values in wave numbers above E-^ 1 

number of outer shell electrons 

|e(i,k) - e(i, j) | in wave numbers 

EZ in electron volts 

B at Te = / (e— -l) 

B at Tr 

He/V 

- LX 

4.803 x 10 , e - elementary charge 

3-l4l6, n 

6.6256 x 10 Planck constant 

10 

2.9979 x 10 velocity of light cm/sec 
I.3805 x 10 Boltzmann constant 

9.1091 x 10 electron rest mass 

2.179 x 10 , Rydberg constant x Planck constant x velocity of light 

8.797 x 10 ^ , ttA , A - Bohr radius 
. 0 0 

number of rows and columns in A matrix 

element of A matrix (in this case I and j refer to serial number 
of level: i.e., A(l5, 23) is ion 1, level 15, and ion 2, level 4 


storage arrays for the indeces of f 


ijk 


12 
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SECTION; 
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SIMULTANEOUS LINEAR EQUATION 
FUNCTION SUBPROGRAM 

NSIMEQ(N,LN,LM,A,B,D,E) 
XSIMEQ(N ,LN,LM,A,B,D,E) 
NDETRM(N,LN,A,D) 
XDETRM(N t LN,A,D) 


TABLE OP CONTENTS 
M or XolMEQ 

PURPOSE 

CALLING sequences 
MEIIIOD 
RESTRICTIONS 
ERROR RETURN 
STORAGE 

miBI M-or XPETRM 
PURPOSE 

CALLING SEQUENCES 

/ 

METHOD , 

RESTRICTION 

.STORAGE 

SUBROUTINE LISTINGS 
TEST PROGRAM 
INPUT DATA 
TEST RESULTS 


M 

1 

2 

2 

J 

3 

4 
4 

6 

6 

6 

7 

7 

8 
11 
U 
IP 




secT,ow = 
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MATH ROUTINES 

1 - - 


UP-3984 .1C 


N5XMEQ NDETRM 122 

XSIMEQ XDETRM 


PURPOSE : 

To solve the matrix equation A*X - B for the unknown matrix X* At the same time, this subroutine 
computes a scaled version of the determinant of the matrix A. 

CALLING SEQVEHCE: 

This subroutine is used with the following calling sequences: 

M = UDIMEQ(N t LN , LM, A f B,D,E) 
or 

M = XSiMHQ(N,LK,LM,A,B,D,E) 

where 

h = a fixed point constant or variable whose value must be equal to the maximum value that 
may be assumed by the subscript I of the matrix A(I,J). This value is identical with 
the value- given in the DIMENSION statement that sets the upper limits for the subscripts 
of A (see sample program). 

LN 2 = a fixed point constant or variable equal to the actual number of rows or columns in matrix 
A and the number of rows in matrix B. 

!J4 = a fixed point constant or variable equal to the actual number of columns in matrix R. 

A - the source program floating point variable used to designate the elements of matrix A, 

• » 

which may be stored by either row or column. The READ and WRITE statements must be 
written accordingly. 

B = the source program floating point variable used to designate the elements of matrix B. 

D = a floating point variable whose value serves as a scale factor by which NSIMEQ or XSIMEQ 

multiplies the value of the determinant of the matrix A. After the execution of the 
subroutine* D contains the scaled value of the determinant. The determinant is formed 
In th» Ini Inwing manner! 

t«N 

’ D » FT A(l»1 )-U 

■ f*t 
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UP-398*. 10 « 

MATH ROUTINES 


SCCTtOM: 

PXC6: 3 


iJSIMLi'v NDEiRM 
Y~.m-U XDETRM 


i.e.. by forming products of success i ve pivots with the proper sign adjustment to com- 
pensate for the row and column interchanges. Unless the value of the determinant is 
needed, it is best to set D equal to zero. 

E — a fixed or floating point variable array of 

length at least equal to the number of rows or columns of matrix A. In this area 

IJ3IMEQ or X3IMEQ keep a record of the column permutations. The contents of E may be 
erased after return from the sub- pro gram. 

M = a fixed point variable which will be assigned the fixed point constants: 

1. if the solution was successful 

2. if underflow or overflow occurs 

3. if the matrix A is singular. 

Cote: T h'; type statement JNTECEE XdJMEQ must appear in the ca L 1 Ing program. 

miu* 

1# Mat be mat icai 

Solution ot the matrix equation A*X » H Is accomplished by upper tr ianqu lar i zing the A 
matrix JSing & maximum pivot for each reduction step. This entai Is searching the reduced 
(N-K+1 )X(!KK+1 ) A matrix - at the Kth stage of reduction - for the element with the largest 
absolute value. A row and column interchange is then performed to bring this element into 
the position. After completion of the triangular ization the X matrix is obtained by 
back substitution. 

2. Coding 

FORTRAN IV for the UNI VAC0 1 07. 

mxmLdm+j 

1« Hie magnitude of N is restricted only by the numh/'r M con’ l'n .itUn 1 . fivuil.iMe and depends 
therefore on the program in which NtllMliQ or XUiMhQ \u used, 

2. U31MEQ and X5IME0 call another subprogram 5IMEQ. Therefore, the use of a subprogram or a 
function named 3IMEQ is forbidden in programs that use NS1MEQ, NDETRM, XSIMEQ, or XDETRM. 

NOTE: To solve the matrix equation A*X := B for the unknown matrix X, the dimensions of the various 
matrices must be: 


A: N x U 

ft: IJ x LM 

X: U x IM 


I.M i N 






iUNIVAC MO"7 BEEF 

MATH ROUTINES 


UP-3984.10 


N3IMHQ :;D£tav 

xsiwe ; joetbm 


It is oossible fcr the programs to give a solution to the matrix equation A * X = B with 
IM > N- if a solution to the matrix equation must be obtained with the B matrix having 
more columns than the order of matrix A* then the dimensions of the various matrices must 


be as follows: 


A: N x LM 

B: N x LM 

x: N x u 



Errors detected in the execution of N5IMEQ or XSIMEQ result in normal UNIVAC Monitor error returns* 

The status of the overflow trigger is not preserved on entry. On exit, this trigger is in the OFF 
position* 

SIQBASt: 

The elements of the matrices A and B are stored in FORTRAN IV order. Execution of this subroutine 
destroys the original A and B matrices. 

After a successful exit from this subroutine, the answers or the X matrix replace the original A 
matrix. This replacement is done according to the scheme: 

A(1,J) Is replaced by X(l,J) 

The programmer may choose to use the same floating point subscripted variable to reference the 
answer that he used to reference the elements of the A matrix, or he may choose to reference the 
solution by a different subscripted variable. If the latter is the case, this variable will appear 
in a dimension and equivalence statement. The choice of parameters which appear with this variable 
in a dimension statement must be chosen in accordance with the following rule: 

If A(I,J) is used to reference elements of the original A matrix and X(I,J) is chosen 
to reference elements of the solution matrix, then' those two variables must appear in 
a dimension statement with the following parameters: 

EQUIVALENCE (A,X) 

DIMENSION A(N,N),X(N,i.) 

where 

J. i number of columns in matrix* 






UNI VAC 1107 BEEF 




UP-3904.10 - * 

MATH ROUTINES 


SECT ION] 

PAGE! b 


NSIMEQ NDETRM 125 

XSIMEQ XDETRM 

STORAGE : 


NSIMEQ requires 41g words of instruction end 4g words of data} SIMEQ requires 740g words of instruc- 
tion and 63 words of data. 
o 

No COMMON storage is used. 
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MATH ROUTINES 

UP- 3984. 10 


NSIMEQ NDETRM 
XSIMEQ X DEISM 


PURPOSE : 

To evaluate the expression 
0 x D£T { A ) 

for a given arbitrary square matrix A and some floating point variable D. 

NDETRM or XDETRM should be used with the following calling sequences: 

M = NDETRM(N,LN ,A t D) 
or 

M = XDETRM(N,LN,A,D) 

where N, LN, and A are as described under NSIMEQ or XSIMEQ and 

D = a floating point variable# Its value will be used to multiply the determinant in accord- 
ance with the scheme outlined under METHOD. On exit from the subroutine the original 
( value of D will be replaced by the scaled value of the determinant. In the case of a 

singular matrix the determinant (D) will be set to zero. 

M - a fixed point variable which will be assigned the following values: 

1, if the subroutine was successful; 

2 * if overflow or underflow occurred; 

3. if matrix A was singular. 


NOTE: The type statement INTEGER XDETRM must appear in the calling program. 

METHOD : 

NDETRM and XDETRM are function subprograms that cause entry into the GiMEQ subroutine. Entry 
by NDETRM or XDETRM cause only the triangular izallon portion of 0JMEQ to be used, (hoe 
NSIMEQ or XSIMEQ for a description of this process). The marmot- in which the determinant 1c 
computed In as follows: During each level of reduction of lire A matrix the expression 


LN 


1=1 


A(X t l) 'Q 
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UP-3984.10 , * 

MATH ROUTINES 


SECTION; 



MSIMC'J iHlF.TRM 127 

XSIMLQ XUFIRM 


io computed. After the final level Of reduction the result iu multiplied by the niv-t, 
Whenever a row or column interchange is required to put a maximum pivot on the diagonal + i i <-• 
sign of the above expression is changed. The final expression will then have the ' e-re. t 
sign* If the second argument of NDETRA/S or XDETRM is a 1 , the error .return M • • w-ll 
indicating the matrix A to be aingular whether or not A ( T , I )~0 . 

2. gfldina 

FORTRAN IV for the UNI VAC 1107. 

R ^TRicnoia i 

1. The magnitude of U is restricted only by t.lv- number of tore locations available and dt-go 1 .'! 
therefore ori the program in which NOJMHQ f >i XCIMliO i * need. 

2. !;C3MisO and XOIMEQ call another subprogram CIMEQ. ITior el oi c , the use •>! a subprogram ur « 
t uric t Ion named U1MEQ is forbidden in programs that use NMMLQ, NIJbTRM, XUIMldj, iV r XI)L.[RM. 

STWCE: 

The elements of matrix A must be stored through the use of a FORTRAN doubly subscripted float inn 
point variable. The subscripts must assume only consecutive integer values. 

NDETRM requires 32 p words of instruction and words of data. 


No COMMON storage is used. 
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PROGRAM FOR SOLUTION OF LINE TRANSPORT 
EQUATION 


i* c 

2* C 

3* c 

1 I* 

5* 

6 * 

7* 

A* 


BOUND * BOUND CASE 

THPEF LEVELS AND CONTINUUM 

SOLVES FOR N3/N1 

DoURLt PRECISION R*l,EMll«G!2«Gt3 

DOUBLE PRECISION R2 , RK, BJ, C12.C21 t C13» C31 ♦ C23, C32» PK i # PK2»PK3, 
lPiK.P2K»P3K»Pl2*P2l,pi5,p3i f p23,p3 2 ,Dl.D2.DJ,PNl?,pH?l ,PH13,PH31 

?PH23»PH3?»ZX*ZY»Eps,ps f Cjt SJ»NL»NlM*WKi,tWKLM#SW,FMl2*FM?l 
DOUBLE PRECISION £F 


f 


^RECEDING 


page bunk not filmed 
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9* DOUBLE PRECISION T 

10 * double precision f 

H* DOUBLE PRECISION FX,GX,AX 

i2* double precision EM*v 

13* DOUBLE PRECISION NOSONI tX ♦ FI *FIX 

14* DOUBLE PRECISION El*E2»E3?E4»XK 

t5* double precision z*zt*dz.st*zkl*con*xni.ri 

lfe* C0MM0 N/BLKi/T( 3I) ,El(3i) f E2(il) *E3(3i) ,E«(31) 

17* C0MM0N/BLK2/F( 31 » 31) ,>1(31*31) ,WL< 31*31) *WLM(31* 31) *Xk (31 *31 ) • 

(8* 1X1(31*31) 

19* DIMENSION EPS (31) *BS(3l)*BJ(31) 

20* DIMENSION EF(31'31)«x(31) 

21* DIMENSION EM(3i,3l) f JCC31),V(2) 

22* DIMENSION WKLM(31*31)*TK(31*5i)fSN(31t31) 

23* DIMENSION XN(3t)*FX(3i*31)*FIX(3l*31)*GX(3i),AX<3l)*RX(3l) *CJC3|) 

?U* DIMENSION WKU3W31) ,SJ(3l) 

?5* DIMENSION ZI31) ,Rl(3i) tRK(3l) IZKLC31) *XN1(31) .XN2C31) tXNK(3l) 

26* DIMENSION DA (2) 

?7* DIMENSION R2(31)»XN3{31) 

?«* DIMENSION Pl2(31)*P2l(Sn»P13C3l)*P31(3n#P23(3l)tp32(3n 

29* DIMENSION XEPS(3l)fETA(3l)*XlOTA(3l) 

30* 10 FORMAT (6E12»8) 

31* 12 FORMAT (1216) 

32* 15 FoRMAT(l6»Dl2,8fEl2,fl) 

33* 20 FORMAT(iH0 f l4*2XtlP5E20#7/(7Xt5E20,7)) 

3«* 25 FORMAT (13H0UNIT PRODUCT) 

35* call date ( 9 » da j 

36* LX * 31 

37* Nc * 31 

36* NR » 31 

39* UMAX ■ 31 

aO* DX * 0,5 

<11* PI * 3.14159265 

<12* TSP * 1,1283792 

a3* ENU * 198305*0 

«4* C * 2*997925EM0 

U5* H * 6,fe256E«27 

46 * HK = 1,38054E*16 

«7* REA0(5*12) IR.IEPs 

U8* READ (5«15)N0tZTfXNT 

«9* N? = 6*ND ♦ 1 

50* NZM * NZ - t 

51* NZP * NZ + l 

52* READ (5*30) (Rl(l) *I»i*NZ) 

53* 30 FORMAT (5015*8) 

54* REAO(5.30> (R2(I) *Xsj*KZ) 

55* READ(5*30) (RK(I).I«i*NZ) 

56* IF(lR,fiT,l) READC5.30) (BJ( I) * **1 »NZ) 

57* READ(5*10) XM.TCUT.VEL 

58* RFAD(5*10) TE*xnE*XNU*C1*G2»63 

59* RF AD(5. 1 0) BK 12* BK2 1 , BK 1 3 * BK31 * BK23 * BK32 

60* READ (5*iO)OMl2*OM2l,OMl3*OM31*OM23*OM32 

61* READ (5*10) B8 12*8 B21»8B13»BB31 ,BB23* 8832 

62* READ (5*10)AP21*AP31,AP32* AB1K,AB2K*AB3K 

63* READ (5*10)AL1*AL2* A| 3*BEi.BE2*BE3 

60 * READ (5*lO)OBlK f OB2K,OB3K.DBl*DB2*DB3 

65* ir(IEPS,EQ,l) READ(5,10) ( XEPS( I ) * ET A( I) * X IOT A( I ) ♦ 1*1 . NZ) 

66* A * (2,0*H*XNU**2)*(XNU/(C**2>) 
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V 7 * 

$9* 

70* 

71 * 

72 * 

73* 

74* 

75 * 

76* 

?7* 

78* 

79* 

80 * 

81* 

82 * 

83 * 

84* 

85 * 

86 * 

87* 

88 * 

89* 


90* 

91* 

92* 

93* 

94* 

95* 

96* 

97* 

98* 


99 * 


100* 

101 * 

102 * 

103* 

104* 

105* 

10** 

107* 

108 * 

109* 


no* 

m* 

112 * 


U3* 

u.«* 

115 * 

116 * 


V H7* 
\l 15* 
ni9* 


120 * 
1 ?1* 


122* 


1?'* 


1 24 * 


HKT * (M*XNU)/(H«*TE) 

Sc « A/(EXP(H«T) » 1.0) 

BET * F.XP(*HKt) 

Cl2 * XNE*0M12 
C2.1 « XM?*0H2l 
Ct 3 « XNE*O m 13 
C31 * XNE*0M3| 

C?3 * XNE*0M23 
C32 * XNE*0M32 

PKl * XNE*(*ll * BEt ♦ XNiE*D8i) 

PK2 ■ *NE*(Al2 ♦ BE2 ♦ XNE*DB2) 

PK3 * XNf*< aLJ ♦ BE3 ♦ XNF*0B3) 

Pi K a AB1K ♦ XNE*Q91K 
PgK * AB2K ♦ X NE*0B2K 
p jK * ABjK ♦ XNE*0B3K 
GO TO (31*36) *IR 

31 D 0 32 I « 1*NZ 

P l 2 ( I) * BB12 ♦ C12 
P 2 1 Cl) * 0021 ♦ C21 * AP2i 
P13CI) * BBl 3 * C l 3 
P3l(I) • BB31 ♦ C31 * AP31 
P23(I) * BB?3 ♦ C23 
P32(l) ■ BB32 ♦ C32 * Ap32 

32 CONTINUE 
GO TO 41 

36 Dn 37 I # 1*NZ 

Pl2(I) * C12 ♦ BK 1 2*BJ( I) 

P?1(I) « Ap2l ♦ CHI ♦ BK21*BJ ( I ) 

P 1 3 ( I ) ■ Cl 3 ♦ BM5*BJ(I) 

P3l(I) ■ Ap3l ♦ C31 ♦ BK31*BJ(I) 

P23(I) = C23 ♦ BK23*Bjn) 

P32(I) * C32 ♦ APS? + BK32*BJ(I) 

37 continue 
«i continue 

Ot a PKl ♦ PK2 

02 a PKl + P*3 

D3 * PK2 ♦ PK3 

PH12 » (P1K*PK2)/D1 
PH2t « (P2K*PK1)/Dt 
PM13 « (PlK*PK3)/02 
PH31 * (P3K*PKl)/02 
Ph23 * (P2K*PK3)/D3 
Pn32 * (P3k*PK 2) Z03 
G12 • G1/G2 
Gl3 a Gl/63 

IFUEPS.EQ,1*AND*IR,EQ*1) 60 TO 43 

DO 42 I a 1*NZ 

EM12 » *(P32CI) ♦ PH32) 

EM21 ■ «*(P23(I) ♦ PH23) 

EmI 1 * P21 ( I) ♦ PH21 • Em 21 
RRl * P12CI) * PHU 

ZX * ( C31 ♦ PH31 ♦ P32CI) + PH32 • (EM12*EM2l /EMI 1 ) ) /aP 31 
ZV « (G13/APJl)*(Ct3 ♦ PH13 - EM2l*RRl/E M ll ) 

EPS(I) * ZX - Z* 

BSCI) a (A*ZY)/EPS(1) 

42 CONTINUE 
GO TO 46 

«3 DO 44 I a t*NZ 
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125 * 
1 ? 6 * 
127* 
1 ? 0 * 
12** 
130* 
ni* 
132* 
133* 
134* 
1 35* 
136* 
137* 
13®* 
13 ** 
lao* 
14 l* 
1 <| 2 * 
143* 
144 * 

145 * 
146* 
147* 
146* 
14** 
150* 
151* 
152* 
153* 
154* 
155* 
156* 
157 * 
150 * 
15 ** 
160* 
1 6 1 * 
162 * 
1 63* 
164* 
165* 
166 * 
167* 
160 * 
16 ** 
170* 
171* 
172* 
173* 
174* 
175* 
176* 
177* 
170* 
17 ** 
1 * 0 * 
101 * 
102 * 


Epsm ■ xeps(i) ♦ etacd 

8sCI) * (XEPS( I)*8C ♦ XIoTA(I>)/EP$(I) 

44 continue 

46 continue 

so * (2.0*HK*TE)/XM + v€L**2 
DNUD * (XNU/C)*SQRT(SQ) 

PHIZ * l,0/(DNUD*SORT(PI)I 

CON « (CH*xNU)/(4.0*P!))*PHlZ*eKl3 

00 35 !*l*NZ 

XNl ( I ) * XNT/(1,0 + P1(I) + R2(I) ♦ RK ( I ) ) 

XN2( I) a R1CI)*XN1(I5 
Xn 3(I) * R2(I)*XNl(i) 

XnK(I) . RK(I)*XNUn 1 

ZkL(I) ■ C0N*XN1 ( I) * ( 1 # 0 - G13*R2( I) ) 

35 Continue i 

ZTl * t«OOOOOl*ZT I 

ZL « ALOGIO(ZTL) 

KL * ZL 

u « KL - ND i 

Z(l) • 0,0 
Z (2) 8 10 ,0**LL 
Z f 3 ) * 2 , 0 *Z< 2 ) 

Z (4) ■ 5,0*Z (2) 1 

NDH * J*ND ♦ 1 1 

NpHP * NOH ♦ J j 

ip cno.eq.d go to 45 1 

DO 40 1*5 » N£)H , 

Z(I) * 1 0 *0*Z ( I<*3) l 

40 CONTINUE I 

45 00 50 J»NOHPfNZ l 

JJ 1 NZ . J M 1 

Z(J) « ZT * Z(JJ) I 

so continue , 

st * 0,0 1 

m> « zd) , 

DO 55 I*2fNZ 1 

DZ » Z(I) «• Z(M) l 

T ( I) ■ ST ♦ 0,5*DZ*(ZKLU*1) ♦ ZKL(I)) \ 

ST * T C 13 I 

55 CONTINUE , 

WRITE (6,60) CT(J) tJ«i,NZ) f 

60 EORMAT(4HOTAU/(1PD3O.10)) I 

WRITE (6,65) > 

65 EoPmAT(19H1BOUnO • BoUNO CA$E* 10 X? 26 HTHREE LEVELS and CONTINUUM , 10 
tX,10HSOLUTION for NJ/N1 ) 

WRITE (6f 70)DA,ND t ZT # XNT t 

70 PORHATdlHOlNPUT DATA* IOOX , A6, A3//20H NUMBER OF DEc*DES *I 3 / 10 H GEf 
lOMETRlC DEPTH *1PE0 ,i/ 10H TOTAL PARTICLES *Eit, 4 ) f 

WPlTE<6*75> XM,vEL 
75 F0RHAT(4H M *|PEll,4/aH V *E11,4) 

WRITE (6, 60) T CUT t 

80 F0RmAT(7H TOUT *F6,1) , 

WRITE (6,05) (Z(I) *R1 (I) t R2(I) *RK(I) tXNKj) ,XN2(I) *VN3(I) *XNK( 

11. NZ) ^ , 

85 FOPMAT(lHO,0X f lHZ t 13X,5HN2/Nl,l0X,5HN3/Nl,10X,5HNK/Nl,ilX,2HNl,13Xf 

t*2HN2,l3X*2HN3,l3X,2HNK//(lP8El5.4)) f 

WRITE (6* 86) (Z(I),ZKL(n.T(I).BJ(I)*EPS(X)*Bs(I) ,I*1,NZ) 

86 FORMAT (1 W0,8X,1 HZ, IPX # 5HKAPPA»11X*3HTAU*12X»4HJBAR,UX»3HEPS*13X,2 
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166* 
1*7* 
1*8* 
18** 
190* 
J 91 * 
192* 
193* 
194* 
195* 
196* 
197* 
198* 
199* 
200 * 
201 * 
202* 
203* 
204* 
2 05* 
206* 
207* 
208 * 
209* 
210* 
211 * 
212 * 
213* 
214* 
215* 
216* 

2 1 7* 

2 1 8 * 

2 1 9* 
220* 
221* 
222 * 
223* 
224* 

225* 

226* 

227* 

228* 

229* 

230* 

231* 

232* 

233* 



237* 

238* 

239* 

240* 


1 H0S/( 1P6E 15.4) J 
WRITE f6.95)TE,XN| f X nu»61*G3 

* E22 - 4 ' 7H ^ 3 ‘ •««.».*- time*,™ *v 

100 0MEG4(1J , j 

,e»s../i*h o«eg»<«J) 

, n _ < 6 * lf, 5)BBl2»BB2!»BBl3»BB3i» BB23 iBB32 ! 

105 FqRMAt(11H BBAr( 12) slpEl6,4/lfH BbarC21) s Fi 6 4/ilH rraomv 

t WRTTE l r6 B |fo)*oRl* E D 6 * <,/1 i H B»i»C23)%E|6.;/uS # Kii?s5? M* J 4 ?' 

lin FS^ATJoi l !2raf? lMP5t ! APi? ’ A01Kf AB2K»AB3K tlt 1 , 

1 10. FqRBAT( 9H Ap (21 ) *1PE18,4/9H 4p(31) *El8,4/9H Apf3»l «FIR ad 

«««*> -ti6.a/ ,h abIrJjJ, Jnf 

us s i rfni i ^^i;nf^?i:j^^PHA, 2 n L.Zl j 

pet * t51 ’ sE,, -“>; 

* T , M D Sn T i£» eu ■«*•«'" »« 

WrITE(M 125) A»BET,8c 

125 F()RhAT( 4H0A s 1 pE23 , 7/4H B sE23.7/9h pUAnCK sfiB 71 

^,10? PE *' EE ^ f EK3»PlK#P2K*P3K*Mt02fD3*PH[2tPH2l,PH13tPH31 f 

lYkVX 

;^m # ?ss h .ks:.j“ , - 4/th phi3 * E2o,<,/7H pH31 

JRITE(6.»») <ZtI).Pl2(I).P21(I),P13(I),P31(I),P23(T).P32n).I«l.NZ 

NT * NZ 


NTM » nT • 1 

9FAO (8*10)YA,yB 

IF (T(MZ) #LE*1 .0) GO TO 440 

XA * DLOGCT(NZ)} 

XmAx e SQRT(XA) 

JX * XMAX/DX 
NX » JX ♦ 2 
XW(l) ■ 0*0 
DO 130 J«2»NX 
XJ * J - l 
XW(J> • 0X*XJ 
130 CONTINUE 

WRITE (6,l35)(XN(KJ,K*ltNX) 
135 FORN AT (5H0X(K) //( IPEj 5*2) ) 
00 140 Isl « NX 
DO HO Jel.NX 
FX(IfJ) * 0.0 
FIX(I»J) * o.o 
HO CONTINUE 

00 145 Is 1 . NX 
FX(Jfl) » 1,0 
Fixn.j) s 1,0 





134 

241* 

145 

i CONTINUE 

242* 


Oq 150 Ja2,NX 

243* 


JM « J * 1 

244* 


00 1 SO 1*1, JM 

245* 


XX * XNCI)/XN(J) 

246* 


FXCI,J) s (1.0 * XX)*(1 ,0 - Y A* XX) 

247* 


FTX(I#J) * FX ( I , J) 

24®* 

150 

continue 

249* 


ie * nosonhfix,x,mx,lx) 

2?0* 


IF (IE.Eq.O) Go To 4?5 

25!* 


WRITE (6,155) 

252* 

155 

Format (Shof (X) j 

253* 


00 160 la 1 » NX 

25«* 


WRITE (6,20) X* (FX(I, j) , Jal ,NX) 

255* 

160 

Continue 

256* 


WRITE (6,165) 

257* 

165 

format ( 1 ihofx inverse) 

25 a* 


00 170 1*1, NX 

259* 


WRITE (6,20) If (FIX( I, J), J*1,NX) 

260* 

170 

CONTINUE 

26!* 


00 175 1 = 1 , NX 

262* 


DO 175 J=1,NX 

263* 


EF(IfJ) = 0,0 

264* 

175 

CONTINUE 

265* 


Do 180 1*1, NX 

266* 


DO 180 J* 1 * NX 

267* 


EF(IiJ) « 0,0 

26®* 


DO 180 Ka 1 » NX 

269* 


EF(IfJ) « EF( I , J) + FX(IfK)*FlX(Kf J) 

270* 

180 

Continue 

271* 


write (6,25) 

272* 


00 185 1*1, NX 

273* 


WRITE (6,20)I,(EF(I,j),J=i tNX ) 

274* 

185 

CONTINUE 

275* 


Gx ( 1 ) a XN(NX) 

276* 


Do 190 J«2 t NX 

277* 


GX(J) * 0.5*XN(J)*(!.o • YA/3,o) 

278* 

190 

CONTINUE 

279* 


WRITE (6,19i)YAfYB 

280* 

195 

FQRMAT(5H0YA «lREu,i,y5H YB «EU,4) 

281* 


DO 200 K*1,NX 

282* 


AX{K) * OtO 

283* 


DO ; 200 Ja 1 , NX 

284* 


A X ( K ) « AX(K) ♦ GX{ J)*FIX( J,K) 

285* 

200 

Continue 

286* 


WRITE (6,205)(XN(I),c,X(I),AX(I),I*t t NX) 

287* 

205 

FoRMATUHOf UX,lMX,l8X»4HG(X)f 17X»2HAK/(1P3E20.8)) 

286* 


NjK * NZ 

289* 


NT KM * NTK m j 

290* 


DO 210 Mai, NX 

291* 


X2 * XN(M)**2 

292* 


Ex2 * EXP (»X2) 

293* 


DO 210 LaltNZ 

294* 


TK(U»M) a T (L) *EX2 

295* 

210 

CONTINUE 

296* 


DO 215 Iat,NX 

297* 


EX » XN ( I ) **2 

298* 


RX(I) * AX(I)*FXP(-Ex) 
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^99* 

215 

CONTINUE 

-Joe* 


00 220 1*1 ,NZ 

301* 


IT * I 

30** 


IF <T(I)»t # 0) 220 f 225*225 

303* 

220 

CONTINUE 

30«* 


ITM * IT 

3$5* 


CO TO 230 

306* 

225 

CONTINUE 

307* 


IT m « IT * t 

3o«* 

230 

CONTINUE 

30<»* 


00 235 I* 1 * NZ 

310* 


is • i 

311* 


IF (T(I)-0,5) 235,235*240 

3 12* 

235 

Continue 

3 1.3* 


ism * is 

319* 


Go TO 245 

315* 

240 

Continue 

316* 


ISM ■ IS * l 

317* 

245 

CONTINUE 

3l8* 


KN * 1 

3{9* . 


DO 250 IM.IMAX 

3?0* 


00 250 JaltLMAx 

321* 


SWtlf J) a 0,0 

3??* 

250 

CONTINUE 

3?3* 

C 

OUTER LOOP ON FREQUENCY — INDEX K 

3?a* 


DO 320 Kal.NX 

32?* 


KK b K 

3g6* 


DO 255 I> 1 « NZ 

327* 


T C 13 * TK(I»K) 

328* 

255 

CONTINUE 

329* 


NTK * NZ 

330* 


NtKM * NtK m 1 

331* 


TsPR * T$P*RX(K) 

332* 


WRITE (6«260)XN(K) tTSPR 

333* 

260 

FORMAT (4HtX *F6,2f iOXfTHCONST eiPEll,4) 

33«* 


WRIT! (6,265) (TK(I,K),I*ltNTK) 

335* 

265 

FORMAT! 1H0.UX»2HTK//(1PE20*7>) 

336* 


Call MMAytNZtTCuT) 

337* 


DO 270 Isl,NTK 

338* 


DO 270 J*1»NTk 

339* 


WKLM(1#J) s WLM(I*J) 

340* 


WKL(IfJ) * NL(I.J) 

341* 

270 

Continue 

342* 

275 

continue 

343* 


IF (IT.EQ.l) GO TO 285 

344* 


DO 280 1st, ITM 

345* 


DO 280 J* 1 , NZ 

346* 


S w C x « J? • SW(I,J) 4 TSPR*WKL(I,J) 

347* 

280 

CONTINUE 

348* 

285 

IF (IT.Eo.nZ) Go To 295 

349* 


DO 290 Is IT , NZ 

\330* 


DO 290 Jsl,NZ 

*s?i* 


SW(I«J) « AMCltJl ♦ TSPR*WKLM(I, J) 

352 * 

‘ 290 

CONTINUE 

353 * 

295 

CONTINUE 

354 * 


WRITE (6,300) 

355 * 

300 

FORMAT(17HOWIJK(L*MBDA « 1)) 

356* 


DO 305 Isl f NZ 



156 


357* 

358* 

359* 

360* 

361* 

362* 

363* 

364* 

365* 

366* 

367* 

368* 

369* 

370* 

371* 

372* 

373* 

374* 

3f5* 

3?8* 

377* 

37»* 

37** 

360* 

361* 

362* 

363* 

364* 

385* 

386* 

367* 

366* 

369* 

3*0* 

391* 

392* 

393* 

394* 

395* 

396* 

397* 

398* 

399* 

400* 

401* 

4 92* 

4.03* 

404 * 

a o 5* 

4q6* 

407* 

408 * 

409 * 

4J0* 

411* 

UJ2* 

4J3* 
4 1 4* 


305 


315 

320 

325 


330 


WRIT E{ 6, 20) X» (WKLM( if J) • J*1 »N7) 

continue 

WRITE (6,310) 

310 FORMAT(13HOWlJK(UMBOA>> 

On 315 Isl'MZ 

WRITE (6 f 20)It(W8LcI f J),J«l,NZ) 

CONTINUE 

continue 
WRITE (6,325) 

FORMAT(7H03(I* J)> 

Do 330 1*1. NZ 

WRITE (6,20)lf (9N(I,J),J*1.NZ) 

CONTINUE 

IF (IT.E0»1) C 0 TO 300 
DO 335 I* | • ITH 

Em<”V^-«».0/U.0 ♦ EPS(I)))*SJ<I.J) 

IF(I.ER.J) EM(I.J) * EMdtJ) ♦ 1.0 
CONTINUE 

IF (IT.Eq.nZ) Go to 350 
DO 345 X*IT*NZ 
00 345 Jal.NZ 

EM(ItJ) ■ •(1.0/EPS(I))9SW(X*J> 4 

IF(I.EO,J) EH(I,J) ■ EN(I*J) ♦ 1.0 '• • ~ ' 

CONTINUE 

Continue 

IF (IT.EQ.l) GO TO 360 
nn 355 TaltlTM 

CJ(I) * (Ep8<I)/U.O ♦ E PS (I)))*BSCI) 

CONTINUE 

IF (IT.EO.NZ) GO TO 370 
DO 365 I*IT.NZ 
CJ(I) * B3(I) 

CONTINUE 
CONTINUE 
DO 37S I*1.NZ 
EMU.NZP) * CJCI) 

CONTINUE 
WRITE (6,380) 

F0RNAT(7H0M(I, J) > 

DO 365 Is 1 * NZ 

WRITE (6t20)Xf(EM(ltj) * J«1 ,NZP) 

CONTINUE 
N a NZ 
MC ■ NZP 
V ( l ) m 4 

CALL GJR(EM»NC,NRfN,HC»84tO, JC,V) 

00 190 1*1 *NZ 

BJ(X) a 8J(X)*a*0 ♦ EP8( I) ) * EPS (I)*BS(I) 

ZY « EPS(I)*B8(I)/A 

It a EPS(I) ♦ ZY 
RRl a P12(I) ♦ PHIP 
EM 12 * -(P32CI) ♦ PH32) 

R?( I )** P (C3 J®!) * (BJ( I) /A ♦ ZY) /( 1 ,0 + BJ(I)/A ♦ ZX) 
JJcn l ( (?ii l Rt ( I) *P2X^*^R2( I ) *P3K) /(PKl ♦ PK2 ♦ 


335 

340 


345 

350 


355 

360 

365 

370 


375 

380 


385 


PK3) 
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0 4 1 5* 
a i6* 
a{7* 
aia* 
a\9* 

4gOf 
4?t* 

4?2* 

4?3* 

424* 

4^5* 

a? 6 * 

«?7* 

4?8* 

4?9* 

«50* 

431* 

4 32* 

433* 

434* 

435 * 

456* 

«37* 

438* 

439* 
aao* 

lAC^OSTTCS 

» 8,23 CPU SECONDS 


SUBROUTINE WMAT(NZ,TCUT) 

CALCULATES W( L AhBDA) and W ( L AmBDA - 1) AT LINE center 
DOUBLE PRECISION WLfWLM 

DOUBLE PRECISION T , E 1 , E2, E3» E4 , A 1 ♦ A2, A3 ♦ A4, AS* ESB, R 1 , B2, B5 
DOUBLE PRECISION F,FJ,XK 

C0MM0N/BLKl/T(3t)»El(31) *£2(31) tE3(3l) *EA(3t) 
C0WM0N/9LK2/F(31#3n »FIC31»31) ,WL{31f3i) ♦WLM(31»3n »XK(31»31J ♦ 
1X1(31*31) 

10 FoPMAT(lHO.l4,2Xf lF5E20,7/(7Xt5E20,7)) 

NT * NZ 
DO 15 I* 1 *NZ 
DO 15 J«1*NZ 
WL(IrJ) * 0.0 
WLM(IiJ) ■ 0,0 

15 continue 

NM « N£ • 1 

DO 25 l» 1 , NH 

JP « I ♦ l 

DO 20 J»JP,NM 

A t * DABS (T ( J-l )*? (!) ) 


t* 

2* C 

3* 

a* 

5* 

6 * 

7* 

8* 

9* 

10 * 

U* 

12 * 

13* 

14* 

15* 

IN* 

17* 

18* 

,9* 

?0* 


390 CONTINUE 

PUNCH 395» (Rl(j) »Jel,NZ) 

PUNCH 395«(R2(n»I«l,NZ) 

PUNCH 395* (RK(J) *J*1,NZ) 

PUNCH 395* (BJ(j) * J» 1 , nZ) 

395 FoPHAT (5F15.S) 

WRITE (6,400)(Z(J)*SJCJ)»BJCJ),R1(J)*R2(J).RK(J),J*1,NZ) 

40 0 FoRMAT(lH0,|3x.tHZ,l9X»lHSf i9X.lHJ t l7X,5HN2/ NJ l,)5X,5H N 3/Nlf 15X*5 Hn 
lK/Nl//(lP6E20,a)) 

ao5 call exit 

410 WRITE (6,415) 

415 FoRmAT(21H0ERR O R i n eo , solving) 

write (6 f 420)NC,NR t N f MC,V 

420 FORMAT (5hONC * 1 3 ♦ 5x » 4HNR »|3*5X*3HN *I3»5X»4 hMC *I3*5X»4hV1 *1^8,1 
1 ♦ 5X ♦ 4HV2 «EU,4) 

GO TO 405 

425 WRITE ( 6 1 430 ) IE* NX 

430 FoRMAT(l4HOERROR IN F(X)/5H IE =i3 ? 5X»4HnX 8 j 3 ) 

WRITE (6,155) 

00 435 I« 1 , NX 

write (6 ,2o)if(FX(i,j), jsi, N x) 

435 CONTINUE 
Go TO 405 

44o WRITE (6,445)T(NT) 

445 FORMAT(20H0TAU-MAX TOO SMALL *lPElt,4) 

GO TO 405 
END 
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21 * 

22* 

?3* 

2 «* 

?S* 

26* 

V* 

28 * 

29* 

30* 

31* 
32* 
33* 
30* 
35* 
3 A* 
37* 
38* 
3^* 

ao* 

at* 

a2* 

«3* 

aa* 

a5* 

a6* 

07* 

08 * 

a** 

50* 

51* 

52* 

53* 

50* 

55* 

56* 

57* 

58* 

59* 

60 * 

61 * 

62 * 

63* 

60 * 

65* 

66 * 

67,* 

68 * 

6** 

70* 

71* 

72* 

73* 

70* 

75* 

76* 

77 * 

7 ** 


IF (M.GT.TCUT) go Tp 25 
A? « DABS(TCJ) * T Cl j 3 
A3 ■ DaB 5(T ( J+l) • T ( 1 ) ) 

Aa » T(J) * T(j-n 
AS a T ( J+ 1 ) * T(J> 

B1 * CT(I)*(E88(At*2) * ESB<A2,2)) ♦ DEXP(*Al) * DFXP<-A?) 

1 ♦ E$B ( A2» 3) • ESB(Ai*3))/Aa 

82 • (T(J*t)/Aa)*(ES8(Alf2) - ESB(A2*2)J - (TC J*l) /AS) * (F.S8C A2* 2) 
1 * CSB(A3f2)) 

B3 a (T(I)*(ESB{A 2,25 - EsB(A3,2)) + DExp(*A?) p DfXP(-A3) 

1 ♦ €SB(A3*3) • E8B(A2,3))/A5 
HL(ItJ) * 0.5*(8l • p2 » B3) 

WlM(IrJ) * WL ( 1 1 J) 

20 CONTINUE 
25 CONTINUE 

Do 35 J«2* NM 

IP * J * 1 

00 30 ialP»NZ 

A3 * OABs'CTd) * T( J*1 ) ) 

IF (A3*GT*TCUT) Gn TO 35 
At * DA8SCTU) * T(J-l)) 

A2 ■ DABs(T(I) • T ( J) ) 

Aa * TCJ) - T(J.l) 

A 5 * T(J+n * T(J) 

Bt * U1/A«)*(ESB(A2,2) • ESB(A1»2)) ♦ C(T(J+1) * T U))/A5)*( 

1 ESBCA3»2) - ESBCA2,?)) 

82 * (0EXP("A3) • DExP A?) + ESBCA2,3) * ESB( A3*3) ) /A5 

B3 a CnFXP(-A2) * DEXPC-'AI) ♦ £SB(Ai,3) - ESP ( A?, 3) ) /AO 

WLdfJ) a 0,S*(B1 + R2 p B3) 

WL m C I* J) » WU(ItJ) 

30 Continue 

35 CONTINUE 

00 00 I»2 * NM 

a 1 * ten • T(j»i) 

A? a T<I*n * TCI) 

B 1 * 2,0 -t ESB( Al f 2) * EsB(A2 t 2) 

82 a CO «5 * OEXPt^Al) * ESB( At *3) ) /Al 

83 a <0,5 ■> DEXPC-A2) ♦ ESB( A2,3) ) /A2 

Wt(IiI) * 0 ,5* (Bl p P2 - 831 

* NU<I»I) • 1,0 

ao continue 
00 85 i*i»nm 
Al a DaBs(T(Nm)*T(I)J 
IF (A1,GT,TCUT) GO Tn 05 
A2 a DABSCTCNZ) - T( I) ) 

AJ a T(NZ) p T(NM) 

Bl a ( T ( I ) * ( ESB( Al , 2) - ESB(A2,2)> ♦ DEXP(-Al) - DfXP(-A2) 

1 ♦ ESB < A2 • 3) - ESB(Ai,3))/A3 
B? a (T(NM)/A3)*(ESB(A1»2) p ESB(A2,2)) 

WLtlfNZ) a Bt . 62 
WlM(I*NZ) a WLCltNZ) 
as continue 

A 3 * T (2) - T(!) 

DO 50 ia2,NZ 
A 2 a 0ABSCTC2) - TCI)) 

IF (A?,6T,TCUT) GO TO 50 
At a OABSCTCl) • TCI)) 

81 « ( ( T ( 2) * T(I))/A3)*(ESB(A?,2) - ESBCA1,?)) 
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79 * 

tti: 
82* 
83* 
8fl* 
85* 
86 * 
87* 
88 * 
89* 
90* 
91* 
92* 
93* 
94* 
95* 
96* 
97* 
98* 
99* 
JO0* 
lot* 
102 * 
103* 
104* 
105* 
1 0 6 * 


02 ■ (D EXP(»*2) • DE*P(*A1) + ESB(A1,3) .. ESB( A2.3))/A3 
WL(I»1) * 0*5*(Bt ♦ 82) 

WLM(I*1) * WU(Itl) 

so continue 

IF (A3.CT.TCUT) GO TO 55 

WL C t * 1 3 S 0 .5* ( 1 *0 • E$B(A3,2) • (0*5 - DEXP(-AJ) + EsBC A3, 3) ) /A3> 
GO TO 60 

55 wuCltl) * 0,5 * C0.25/A3) 

60 WIM<1»1) « WL(ift) * 1.0 
A i « T(NZ) * T(NM) 

IF (Al.GT.TCUT) GO TO 65 

WL(NZ*MZ) * 0 »5*( 1,0 • Es»(Al,2) - (0.5 - DtxP(-Al) ♦ EsB(Al»3)) 

1 /*1) 

GO TO 70 

65 WL(NZfNZ) ■ 0,5 » (0,25/Al) 

70 WlM(MZ.NZ) * WL(NZ.NZ) * 1,0 
WRITE (6,75) 

75 FORMATCIOHOW(LAMBDA)) 

Do 80 I=1»NT 

WRITE (6, 10) I. ( WL( I.J) ,J»1 , NT ) 

80 CONTINUE 

WRITE (6,85) 

85 FORMAT(14HOW(LAM8DA - 1)) 

Do 90 I»1 »NT 

WRITE (6,10)I,(WUH(I ( J),J*1,NT) 

90 CONTINUE 

RETURN 

END 


iagnostics 


T!pN TIME * 2,60 CPU SECONDS 


l* 

2 * 

3* 

a* 


5* 

6 * 


7* 

8* 

9* 


10 * 

11* 

12* 

1.3* 

1«* 

15* 

1 ** 

if* 

10 * 

19* 

? 0 * 

21* 

22 * 

23* 

2«* 

25* 


C 

C 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


SUBROUTINE GJR( A»NC»NR,N,HC»$, JC.V) 
DOUBLE PRECISION A,X,V 
DIMENSION A(NR,NC),JCU)t v (25 


.?c is the permutation vector 

KD IS THE OPTION KEY FOR DETERMINANT EVALUATION 
KI IS THE OPTION KEY FOR MATRIX INVERSION 
L IS THE COLUMN CONTROL FOR AX«B 
M IS the COLUMN CONTOL for matrix INVERSION 


INITIALIZATION 


IWPVM) 

Me 1 

s*n 

L«N*(MC*N)*dW/4) 
KD*2*MOD( IW/2,2) 

IF (KD*EO, 1) V (2)*0, 
K I»2*M00( IW » 2) 

Go TO (10»20)»Kl 


initialize JC for invfrsion 


10 DO 15 I*1»N 
15 JCCI)»I 
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26 * 
27* 
28 * 
2^* 
30* 
31* 
32* 
33* 
3«* 
35* 
36* 
37* 
3«* 
59* 
4) ft* 

at* 

a2* 

a3* 

«a* 

a 5* 
a 6* 
a7* 
a6* 
a9* 
50* 

51* 

52* 

53* 

5«* 

55* 

56* 

57* 

58 * 

59* 

60 * 

61 * 

62 * 

63* 

60 * 

65* 

66 * 

67* 

68 * 

69 * 

70* 

71 * 

72* 

73* 

ja* 

75* 

76* 

77* 

78* 

79* 

80* 

81* 

82* 

83* 


C 

C 

C 


e 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


e 

c 

c 

c 


SEARCH for PIVOT row 


25 

30 


35 


«0 


Do 85 1* 1 # N 
GO TO (30*25) *KI 
mb! 

IF (I.ea.N) GO TO 55 
X«M • 

DO 35 jBltN 

IF (X,GT,ABS(A(J»I))) GO TO 35 
X « DAPS(A(J»i)) 

K*J 

CONTINUE 

IF (K.FQ.I) GO TO 55 
S*»*5 

V(I)*»V(I) 

GO To (ao»a5)*KI 
MU«JC(!) 

jc( n«jetK) 

JCCK)*mU 



INTERCHANGE row I AND ROW k 


as DO 50 J*m,l 
X c A ( 1 1 J) 

A(I,J)bA(K, J) 

50 A(K,J)bX 

test for SINGULARITY 
55 IF CDABSCA(I»m «GT.o*0) GO TO 60 


matrix is Singular 


IF(KD.FG.I) V(1)b0, 

JC(1)*I-1 
RETURN 6 

60 GO TO (65» 70) *KD 

Compute the determinant 


65 IF(A(I» I) .LTiO,) s*«s 

V (2) * V ( 2) + DlOG(DABS<AU*I))) 

7ft X*A(I*I) 

An*i)*i. 



REDUCTION OF THE I-TH ROW 


00 75 JbM*L 
A(I*J)»A(I*J)/X 

TEST OVERFLOW SWITCH, IF ON 
RETURN NEGATIVE VALUE OF I IN JCC J ) 

CAUL OVERFL CIFL) 

IF UFL.EQ.l) GO TO 120 
75 Continue 


c 



e«* 

ss* 

06* 

07* 

AS* 

09 * 

90* 

91* 

92* 

93* 

94* 

95* 

96* 

97* 


90* 

99* 

too* 

lot* 
102* 
103* 
10 ** 
10 5* 
!06* 


107* 

10®* 

10 ** 

110 * 

lit* 

112* 

113* 

11** 

U5* 

116* 

117* 

1 ) 0 * 

119* 

120* 

1?1* 

i?2* 


1?3* 

12 «* 

125* 


C REDUCTION OF ALL REMAINING ROMS 


00 85 K*1 • N 
IF (*,EQ,I) GO TO 05 

X*A(K,i) 

A(KtI)*0. 

00 80 J* M *L 
A(K,J) B A<K f J)-X*A(I # j) 

C * 

f TEST overflow SWITCH, IF on 

c RETURN NEGATIVE VALUE OF I IN JC(1) 


CALL OVERFl flFL) 

IF (IFL,CQ,t) GO TO J20 
00 CONTINUE 
85 continue 

c * — ----- 

C AX*0 AND DET, ( A) ARE NOW COMPUTED 


GO TO (90 * 1 lS) #Kl 

C 

c PERMUTATION OF THE COLUMNS FOP MATRIX INVfRSION 


90 D 0 110 Jsl.N 

IF (JC(J).EQiJ) GO TO 110 

JvT*J+l 

D 0 95 I«J4,N 

IF ( JC( I) *EQ.*J) GO TO 100 
95 CONTINUE 
100 JC(I>«JC(J) 

DO 105 K*1 f N 
X*A(K*I) 

A(K,I)*A(Kf J) 

105 A(K»J)*X 
110 CONTINUE 
11$ JCCI)*N 

IFCKD.EQ.I) V(l)*S 
return 

120 JC C 13*1-1 

I F ( Kp, EQ « 1 ) V(1) = S 

RfTURN 6 

EnD 


I AGNOSTICS 


T I ON TIME * 


2,00 CPU SECONDS 



1* 


OOUBU PRECISION FUNCTION 

2* 


OoUBtE PRECISION A,X,F 

3* 


DIMENSION A(1),X(1) 

a* 


N « L m ! 

5* 


max * N*LHAX ♦ L 

9* 


MAX « N*LMAX ♦ l 

7* 


oo to i«ta 

8* 


X(I) * 1,0 

9* 

10 

continue 



K1 » * LMAX 

11* 


00 55 K*1 » L 

12* 


*1 « *1 ♦ LMAx 

1** 


K2 • K l + K 

1«* 


IF CACK2)) 15,80*15 

J5* 

15 

PO 30 1*1, L. 

19* 


Jl • M ♦ I 

17* 


IF (A(J1)) 20,30*20 

|W 

20 

F ■ 1 #0/A( Jl) 

J9* 


X 1 1) a X ( I ) * F 

20* 


DO 25 Jlal, MAX, UMAX 

21* 


A(J1) a A<J1)*F 

22* 

25 

CONTINUE 

23* 

30 

CONTINUE 

29* 


A(K2> a X (K) 

2*5* 


X (K) a l t 0 

29* 


DO 50 Ibi,l 

2** 


HI « K • I 

28* 


IF CKIJ 35,50,35 

29* 

35 

Jl '■ K| * I 

30* 


IF ( A( Jl ) ) 00,50*00 

31* 

00 

A(J1) a 0,0 

32* 


DO 05 J2*I, MAX, UMAX 

33* 


Jl ■ J2 + K I 

30* 


A(J2) a A(J2) - A(Jn 

33* 

05 

continue 

36* 

50 

Continue 

37* 

55 

continue 

36* 


00 70 Ial,N 

39* 


IF (Xtm 60*80*90 

00* 

90 

F a 1,0/xCi) 

ill* 


DO 65 J1»I,MAX,LHAX 

02* 


AfJt) a A ( Jt ) *F 

41 3* 

95 

Continue 

40* 

70 

continue 

415* 


NOSONI a 1 

09* 

75 

RETURN 

07* 

80 

NOSONI * 0 

tie* 


G 0 TO 75 

09* 


end 


NOSONHA, X»l*LMAX) 


JiAGNOSTicS 

'7 ION TIME ■ ,93 CPU SECONDS 
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IMFUT DEFINITIONS 


General 
IR = 

IE PS = 

ND = 

ZT = 

XNT = 

R1 = 

R2 = 

RK = 

BJ = 

m = 

VEL = 
TOUT = 

TE = 

XNE = 
XKU = 

G1 = 

G2 = 

G3 ^ 

YA = 

YB = 


1 at initial start 

2 for iteration 

1, e^, and are read in from P48 (Code l) 
IEPS f 1, T] and are calculated 
number of decades (max. = 5) 
geometric thickness 
total number density 

n 2 / n i 

n 3 / n i 

y-i 

J (input on iteration) 

^ Ph. 

mass of Helium = 6.6408 x 10~ g 

non- thermal velocity in Doppler width 


upper limit on optical depth to avoid overflow in exponential 

routines 

T e (°K) 

n e 8 

v i j(l/sec)( = c x 10 from program 48 (Code l) 


statistical weights for each level 


constant coefficient in / (x) . and g . ( usually = 1. 

tJ 0 

constant coefficient in /' (y) . and g’ . (usually = 1. 

0 J 


0) 

0) 


Two Level Atom (Bound -Free Case) 


DY = 
BIJK = 
BIKJ = 
ENU = 
ALPHA = 


increment in y ^ 

B^ Einstein absorption coefficients - Ion I 
B^ Einstein absorption coefficients - Ion I 
AE^ in wave numbers 

^2 



Two Level Atom (Continued) 


BETA = 

^k2 


ABAR = 

A 2k 


OMEGA = 

Q 21 


OBAR = 

s* 


DEAR = 

«*> 


AP = 

A 21- 


BB = 

B v B 21 


RS = 

hk at 


OB « 

D Ik 


XE1 = 

* 

* 


XE2 = 

n 2 


XEK = 

* 

”k 


0M12 = 

°12 


BB12 = 

B B 
v 21 


Input Card 

Formats 


Card 

Format 

Content 

1 

( 16 ) 

IR 

2 

( l6 , 2E12 . 8 ) 

HD, ZT, XMT 

3 

(6E12.8) 

TCUT 

4,4a, etc. 

(5E15.8) 


5,5a, etc. 

11 

R Ki 

6,6a, etc . 

fl 

BJ^ (only on iteration) 

7 

( 6E12 . 8 ) 

BY, BIJK, BIKJ, EHU 

8 

II 

TE, XNE, ALPHA, BETA, ABAR, OMEGA 

9 

IT . 

OBAR, DEAR, AP, BB, RS , OB 

10 


XE1, XE2, XEK, 0M12, BB12 

11 

1! 

YA, YB 

Three Level Atom (both 

xij n, and nJn. solutions) 



Einstein absorption coefficient 



14-5 

Three Level Atom (Continued) 

BB. . = 

B B, . 



v kj 


AP. . = 

A. . 


1J 



AB. = 
lk 

A ik 


AL. = 

l 

°k,i 


BE. = 

B, • 


i 

_ k ,i 


OB, 1 = 

q.. 


ik 

_ik 


DB. = 

a . 


l 

k,i 


Input Card 

Format 


Card 

Format 

Content 

1 

( 216 ) 

IR, IEPS 

2 

( l6,D12.8,E12.8) HD, ZT, XNT 

3. 3a, etc. 

( 5 D 15 . 8 ) 

Rl. 

1 

4, 4a, etc. 

TT 

R2. 

5,5a , etc. 

Tl 

1 

RK, 

x 

6, 6a, etc. 

1! 

BJ. 

i 

7 

(6 E12.8) 

XM, TOUT, VEL 

8 

r» 

TE, XHE, XNU(l,j ) , Gl, G2, G3 

9 

ir 

BK12, BK21, BKl3, BK3l, BK23, BK32 

10 

n 

0M12, 0M21, QM13, QM31, OM23, OM32 

11 

it 

BB12, BB21, BB13, BB3l, BB23, BB32 

12 

Tl 

AP21, AP31, AP32, ABLE, AB2K, AB3K 

13 

tt 

AL1, AL2, AL3, BE1, BE2, BE3 

14 

Tl 

OB1K, 0B2K, OB3K, DB1, DB2, DB3 

15 

Tl 

XEPS, ETA, IOTA, XEPS, ETA, IOTA 

16 

IT 

YA, YB 



1 U6 


omega 

AP 

EB 

pBAR 

DEAR 

ABAR 

pD 

HH or H 


Q. , electron collisional excitation coefficient 

ljk 7 

A. „ , Einstein spontaneous transition probability 

1JK 

BB. . , Einstein absorption transition probability x Planck function 
13 K 

27 . , electron collisional ionization coefficient 

if. , electron collisional recombination coefficient 
K. radiative ionization rate 
optical depth 

physical thickness of layer (calculated in AMAT) 


13 



l J +7 


DEFINITIONS OP OTHER CONSTANTS AND VARIABLES 


PI = 

C = 

H = 

HK = 

C « 

P 

TSP = 
Z = 

DZ = 
T = 

NT - 
TK = 
NTK = 
X r= 

DX = 
EX = 
FX = 
FIX = 
GX = 
AX = 

Y = 

DY = 
NY = 
FY = 
FIY = 
GY = 
AY = 
XN1 = 
XN2 = 
XN3 = 
XEK = 
WL = 
WLM = 


it = 3.14159265 
2.997925 x 10 10 cm/sec 

6*6256 x 10 erg s y Planck constant 
"■l6 o 

1. 38054- x 10 erg / K, Boltzmann constant 

1.8 X 10" 18 , 7-6 X 10" 18 

1.1283792 = 2/A 

geometric depth 

increment in Z 

table of optical depths at line center 
number of T 

optical depth dependent upon frequency- 
number of TK 

( v " v 0 )/( Av d ) 
increment in X 

number of X 


/ (X). 

J 

inverse matrix of / ( X) 

v Aia 

increment in Y 
number of Y 

inverse of / T (y) „ 

J 

A',. 


n l 

°2 

n 0 


n. 


(A) 


W at line center 

w. A) 

1J 


at line center 



l4S 


WKL = 
WKLM = 

El = i 
E2 = 

E3 = 

E4 = , 


P. . = 
iJ 

PH. . 
ij 

EPS = 
BS = 
EA = 
EB = 


EM 


1J 


W. ^ at frequency v, 

ljk (A-l) k 

at frequency 


exponential integrals of orders 1, 2 , 3 and 4 


refer to capital P^ 

refer to script P. . 

ij 

e 

B S 

a 

0 

b 

e 

M ll ? M 12 ; etc ‘ in ^ erivil| S X and Y in multi-level case 




7.7. OJR 




- Oetoniiinunt; Inver*®; Soh;fi</n of S J truihoncou r . Equations 


’.s: pCSO 

"his subrouti v,r solves simultaneous equations, computes a determinant, or 
diverts a matrix or any combination of the three above by using a Gauss- 
bj r d a n o 1 im ina t ion t o c hn i * j u e vrl th c c 1 urnn pivoting. 


; s „ r Procedure 


Entry 


where: 


DESCRIPTION 


TYPE 


A 


NC 

NR 

N 

MC 


k 


JC 


is the matrix whose inverse or floating-point array; 

determinant is to be determined. input and output 

If simultaneous equations are to 
be solved the last MC - N columns 
of the matrix are the const ant vec- 
tors of the equations to be solved. 

On output if the inverse is computed 
it is stored in the first N columns 
of A. If simultaneous equations are 
solved, the last MC - N columns 
contain the solution vectors. 


is the maximum number of columns of 
the array A. 


FORTRAN integer; 
input 


is the maximum number of rows of 
the array A. 


FORTRAN integer; 
input 


is the number of rows of the array 
A. 


FORTRAN integer; 
input 


is the number of columns of the FORTRAN integer; 

array A* This entry is a dummy input 

argument, if simultaneous equations 
are not solved. 


is a statement number in the calling input 
program to which control is returned 
if an overflow is detected. It must be 
preceded by $ in the calling sequence. 


is a one-dimensional permutation a: ray FOR 1'RAN integer array; 

of N elements used for permuting the input and output 

rows and columns of A if an inverse 

is comoutod. If an inverse is not 

c om pu ted this a r qua: o n t m u s t have at 

legist o fu- cel! ft# r th e ; ' r r o r re t u r n 

identification. On output, the first 

element of the array in N if control is 

re turned uonnally. j r an overflow is 


• U r»-7,S42 


UN I VAC HOB WiATH 


PACK 


SECTKiMt 
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dp tec ted, the first element is 
the negative of the last correctly 
completed row of the reduction. 

If matrix singularity is detected,' 
the entry contains the value of 
the last row before the singularity 
was detected . 

V is a one-dimensional array. If the floating-point array} 

determinant is not computed it has input and output 
one entry, otherwise it has two. On 
input V(1 ) is the option indicator, 
its values are set as follows: 


Operation 

1. 

2. 

V(1) 

3. 4. 

5. 

6. 

7. 

Compute Determinant 

no 

yes 

yes 

no 

no 

yes 

yes 

Invert Matrix 

yes 

no 

yes 

no 

yes 

no 

yes 

Solve Equations 

no 

no 

no 

yes 

yes 

yes 

yes 


On normal return from the program 
V(1) contains the value of the 
natural logarithm of the absolute 
value of the determinant and V(2) 
contains the sign of the deter- 
minant. If an error return is 
made, and the determinant was to 
be computed, then V(1) is set to 
0 and, if an overflow return was 
made, V(2) contains the sign of 
the last correct partially-computed 
value of the determinant. 


7, 7. 2. 2. Restrictions 
None. 


7. 7. 2. 3. Special Considerations 

(1) If the matrix is singular or ill-conditioned, roundoff error may 
cause large discrepancies in the results, 

(2) In the case of a singular matrix, return may not be made through 
the singularity exit because of roundoff error, 

(3) See paragraph 7. 1.2. 3. for notes on usage of the row-dimension in 
arguments N and NR. 


7. 7. 2. a. • Other Subprograms Required 


DACKl 71 


None 
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7. 7. 2. 5* larror Returns 

(1) If a singularity is detected, the first element of the array JC is 
set to the row number before the singularity was detected and, if the 
determinant was to be computed the value of V(1) is set to 0.0. Control 
is then returned to the calling program at the statement number speci- 
fied. 

(2) If an overflow is detected, JC(1)is set to the negative of the last 
correctly completed row of the reduction. V{2) is set to the sign of 
the partial value of the determinant that was computed until this 
time. Control is then returned to the calling program at the statement 
number specified. 


7.7.3. Supporting Information 


7. 7. 3.1. Mathematical Method 

For any matrix A, if a matrix B exists such that BA - AB = I and I is 
the unit matrix, then B - 

If AX = C, where A is n by n, X is n by p, and C is n by p, then the 
solution to these sets of simultaneous equations is 

X = A _1 C 

The determinant of A is defined by the following equation 

A| = V n a . 

i=1 J i 

where 

is the (i,j)th element of the matrix A. 

f( j-j » • • - • j^) is the number of transpositions required to transform 
to J n )> summation is over all permutations 

of the integers (1,... f n). 

The solution to all of these problems is found by using a Gauss- Jordan 
elimination scheme with row scaling and maximal pivoting by columns. 

Faddeev, D.K., and Faddeeva, V.N., Computational Methods of Linear 
Algebra , W.H. Freeman and Co., (1963). 

Ralston and Wil f. Numerical Methods for Digital Computers Wiley (i960). 


7*7.3. 2. Programming Method 

(1) For each column below the diagonal, the program searches for a pivotal 
element by finding the element of maximum absolute value in the remain- 
ing rows of the column. 

« 

(2) This row is interchanged with the row of the diagonal. 


'■-V 

7.7. 

7.7. 


7.7. 

7.7. 


i isshvag moo 
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(3) Each of the elements of the pivotal row is divided by the pivot except 
the pivot which is replaced by its reciprocal. 

(4) All the other rows of the array are changed by the formula 


a . . ~ a. , — a . , a. , 
ij ij lk kj 

where is the pivotal element. If i - k, a^ is replaced by 0. 


(5) When this process has been completed for each diagonal of the array, 
the columns of the matrix are repermuted to give the inverse in the 
first N columns of the array A. 

(6) If the determinant is to be found, each permutation of rows and columns 
changes the value of its sign. The natural logarithm of the absolute 
value of the diagonal element is summed after step 2. 


(7) Only the computations necessary for the options specified are carried 
out. 


3.3. Storage 

GJR » 470 positions, not including the Library Subroutines OVERFL, 

ALOG and NERRS2. 


3.4. Nomenclature 

IFL is the overflow test indicator. 

IW is the FORTRAN integer value of the option indicator. 

KD is the option key for determinant evaluation. 

KI is the option key for matrix inversion. 

L is the ^olumn control for solution of simultaneous equations. 

M is' the column control for matrix inversion. 

MU is a dummy variable. 

S is the sign control for determinant evaluation. 

X is a dummy variable, 

• All other variable names used in GJR are either loop indices or are 
defined in paragraph 7, 7,2.1. 


4. Test Design 


4.1, Introduction 

The test program reads a matrix and seven options into main storage 
and successively references GJR for each of the options, thus testing 
the program's functioning. 





iVAC n n rj :r: ; 


n r 7 


7. 7*4. 3. 


Common l 

The test program computes s 

A" 1 A = R ~ I 

after the inverse ir> found by GJR. This check shows that computational 
accuracy is bettor than A x 10”^ for all elements of the inverse matrix, 
however, this accuracy depends on the well-conditioning of the matrix 
to be inverted. 

Test Input 

The matrix used as a test matrix was symmetric positive definite of 
order A scaled so that each diagonal element is. one. See Fadeeva 
Chapter 2. 


Card 1 : 


Card 2 j 


Title Card, FORMAT (12A6) 

Col, 1-72 An alphanumeric heading to be 

printed by the test program. 

Parameter Card, FORMAT (312,715) 

Col. 1-2 N , the number of rows of matrix A, 

Col. 3-4 NO, the number of columns of matrix 

A, 

Col. 5-6 NOP, the number of options to be 

specified in successive calls to 
GJR. NOP < 7. 

Col. 7-11 I OP ( 1 ) , first option indicator. 

Col. 12-16 \ IOP(2) second option indicator. 


Col. 71-80 


IOP(NOP ) , last option indicator. 


Each element of the I OP array appears as V(1) in a separate call to GJR. 


Card 3f f i 


Matrix Cards, 
Col. 1-10 . 

Col. 11-20 f 


Col. 71-80 ' 

The input data is listed below. 


FORMAT (8E10.3) 

Elements of A. The NC 
columns of row 1 are input 
first, followed by rows 2 
through N. 


TEST MATRIX 1. 

(PEP. FADDtLEVA 

- C0MP.MTHS.0F 

040607 

t 

2 3 4 

5 6 7 

i. 

.42 

.54 

.66 .25 

.42 

1, 

.32 

.44 .45 


,32 

2. 

.2? .65 

,66 

.44 

.22 

1. .8? 


<v & \D cr 
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For inputs, definitions, and subroutines see Appendix B. 


1* 

2 * 

5* 

a* 

*5* 

6 * 

7* 

8* 

9* 

io* 

n* 

12* 
H* 
m* 
15* 
16 * 
1 7* 

J ft* 

t 9 * 
20 * 
21 * 
22 * 
23* 

2«* 

25* 

26 * 

27* 

pfi* 


C. HOUND - FREE CASE 

c two levels amp continuum 

DOUBLE PRECISION FY*FIY 
DOUBLE PRECISION T »£ \ , E2*F3*E4 
DOUBLE PRECISION F*Fl»XK*EFtX 
double precision Em . v 

DOUBLE PRECISION Z » 2y * HZ ♦ ST • Z«L» XN } 

CoHM0N/8LKJ/T(31)«Fi(31).E2C3t)*E3C3l)»E4(3l) 

C 0 MM 0 N/&lK 2 /F( 3 t» 3 l),Fl( 31 f 3 n f WL( 3 l, 3 i),WL.M( 3 i, 3 n ,XK( 3 l, 31 ) f 

1 XL C 3 1 » 3 1) 

DIMENSION RPC3n*EMf3i f 3l).B(3i) *JC(3I)*V12) 

DIMENSION EE(3|»il)tX(3n 

DIMENSION wKLMC3t»31)*TK(3i.3l)tSN(31f 31) 

DIMENSION Y(3l),FYC3t,31),FIYC31»3l).GY(31),AYC3!) 

DIMENSION Z(il)*Rl(3l)*RK(31)»ZKL^31)»XNl(31)»XM2(3l)fXNKC3l) 

DIMENSION DA (2) 

DIMENSION BJC3l)*P12f3n»P2lC3l>»EA<3n*te(3l) 

DIMENSION WP(3J * 3 1 > , WKL (31*31 ) ,RA(3t) 

11 F0RMAT(I6*D12,8»E12,r) 

16 Format (6e12*®) 

21 F0RMAT(I6.2E12*B) 

26 FORMAT ( lHO,Ift, 2 Y»lP< 5 F 20 * 7 /( 7 X t SE 20 , 7 )) 

31. FORMATC13HOUNIT PROD(iCT) 

CALL DATE(9»0A) 

DO 36 1=1 *|,MAx 
DO 36 J=1 »L MAX 
wp(I»J) s 0,0 
SW(IiJ) = 0.0 


P-RJdCEDlNG PAGE BLANK NOT 


FILMED 
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?9* 
30* 
31* 
32* 
33* 
3« I* 
35* 
36* 
37* 
38* 
3 9* 
« 0 * 
al* 
a2* 
a3* 
aa* 
as* 
a&* 
# 7* 
a8* 
a9* 
SO* 
si* 
52 * 

S3* 

5«* 

SS* 

56* 

S7* 

5«* 

59* 

60 * 

hi* 

62 * 

63* 

hH* 

65* 

66 * 

67* 

68 * 

69* 

70* 

71* 

72* 

73* 

79* 

75* 

76* 

77 * 

78* 

79* 

80* 

81* 

8 ?* 

83* 

8a* 

85 * 

86 * 


36 CONTINUE 
NC » 31 
MR * 3l 
l BAX S 31 
C * 2,9979?Sfc*lO 
H * 6.6256E-27 
HK a 1.38059E-16 
PI a 3 * 1 a 1 5926 
CP * 7.6f-t6 
READ (5*21)IR 
READ (5»inNl?,ZT*XNT 
READ (5*16)TCUT 
MZ a 6*ND +1 

NZM S.Kil m l 

NZP = NZ + 1 

RpAD (5. abb) <R1(D *Ist »NZ) 

READ ($*abb> CRK(i) »Iai »NZ) 

IF CIR,GT.U RFAO (5, «bb) ( B J ( I) ♦ I* 1 * NZ ) 

RE AOC5» lb) OY.filJKf BIKJiFNij 

READ (5*1 b)TE t XNE» ALpHA * RfcTA »ABAP»Q mFGA»u» AR *08 ar,ap* BP* RS»08 

READ (5* tb)XEl*XE2*XFK,0Ml2,BBt2 

CON * (8,0*Pl)*( CENU**3)*CP)*C 

ThN a H*ENO*C 

THD a MK*TE 

TH a THN/THD 

C 1 K a XNE*OB 

P? K s ABAR 7 XnE*08Ar 

PK2 a ALPHA * SETA + XNE*DBAR 

PK2 s XNF*PK2 

GO TO caps!)*!* 

a i oo ab i» i * mz 

P 1 2 ( I) a BP 1 2 + XNE*0M12 
P? t ( I) a BB ♦ AP + XnF* OMEGA 

ab continue 
go to 61 
51 r>0 5b I« 1 * nZ 

Pj2(j) a BIKJ*BJ(D ♦ XNE*0M12 
P 2 1 Cl) * BIJK*BJCI) + AP t XNE*OMEGA 

5b continue 

61 on 66 I=1*NZ 

Pm S P?l(I) 7 P2K 
RB a (Pi2U)*P2K)/PN 
PA * (P21 (I)*PK2*(XEK/XE1))/PN 
F. 8 ( I ) a (ClK + P«)/R5 
Ea(I) a (ClK ♦ PA) /Rs 
66 CONTINUE 

OO 71 !*1*NZ 

XN1(I) = XNT / ( 1 *0 + R 1 ( I ) 7 RK ( I ) ) 

XN2( I) a R1(I)*XN1(I) 

XnK(I) a RK ( I ) *XNl C I ) 

ZKL(I) a CP*XN1(I) 

71 CONTINUE 

ZlL a l.OOOOOlfZT 
Zi a ALOGIO(ZTL) 

Kl * ZL 
U. a KL - NO 
Z(l) = 0.0 
Z(?) a 10.0**LL 
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90* 

91* 

9?* 

93* 

pa* 

9S* 

9ft* 

97* 

9ft* 

99* 

1 no* 
1 n l* 
1 o2* 
1 n3* 
1 04* 
10S* 
10** 
1 0?* 
1 oft* 
1 09* 
1 1 0* 

111* 
1 13* 

1 1 3* 

1 1 4* 
1 IS* 
116* 
1 i 7* 
1 1 ft* 
1 1 9* 
1 20* 
121* 
1 ?2* 
!?3* 
1 ?a* 
12S* 
126* 
1 27* 
1 2ft* 

1 ?9* 
1 30* 

131* 
132* 
1 33* 
1 34* 
1 35* 
1 3ft* 
1 37* 



141* 
1 42* 
143* 

1 44* 


2 ( 3 ) s 2 , 0 *Z( 2 ) 

2(4) = 5,0*Z(2) 

MDH = 3 * NO t 1 

NOWp s NOH * 1 
IF (NU.EQ.i) GO TO 
DO 7 ft I*StKDH 
Z(I) = 1 0 , 0 *Z ( 1 - 3 ) 

7 6 CONTINUE 

ftl Do 66 J*NDMP»NZ 

JJ s MZ . J + i 
Z(J) s ZT • Z(JJ) 

66 Continue 
ST * 0,0 
T ( 1 ) = Z(l) 

DO 91 1*2 * NZ 

D7 * Z ( I ) - Z( j-1) 

T ( I ) : ST + O.S*DZ*C7KLCI-l) * ZKL(I)) 

ST » TCI) 

9i Continue 

WRIT f. (6,96) (TCJ) , J»t *NZ) 

96 FOPMaT(4h0TaU/(IPD30.16) ) 

WRITE (6,101) 

ioi porhatobmibuuno * free case*iox» 2 <iht*«o levels and continuum) 

WrItF ( 6 f l06)DA*NDtZT*XNT 

106 FnPMAT(nM 0 INPUT DAT A f 1 0 0* » * 6 , A 3//2 0 M NUMBER OF DECADES 3 I 3 /I 8 H GE 
lOMETRIC DEPTH =lPE6,i/18H TOTAL PARTICLES =£11,4) 

WRlTt(6*lin CP*DY 

111 F0RMAT(5H CP * IPE1 1 »0/5H DY =Eft.l) 

WRITE (6* 1 16)BlJK f eiKJ»ENU 

116 FORMAT ( 7 H BIJK ■lPEH.a/ 7 H 9 IKJ sEU. 4 / 6 H NUK =F 12 . 5 ) 

write ( 6 , 121 )tcut 

121 FoRMAT( 7H TCUT =Fft t l j 

WRITE (6,126) (Z(I)»Rl(I)»RK(I),XNlCl)*XN2Cl),XNK(I),ZKLCl),T(I)f 1 = 
1 1,NZ) 

126 FORMAHlHO, 6 X f 1HZ,13X, 5HN2/N1 , 10X*5HNK/N1, UX,2HN1 , 1 3X . ?HN?, 1 3X » ?H 
1NK,12X,5HKAPPA, ilX,3HTAU//(iP8El5.4)) 

WRITE (ft, 136 )Tfc,XNE,nMi 2 ,OMEGA,OB*OBAR,DeAR,RSf ABA r,BH| 2 , 0 S*AP, ALP 
1 HA» 8 ETA 
WRITE (6,131)TH 

131 FORMAT (Bh THfcTA *1PE19,4) 

WRlTE( 6 , 132) CON 

132 FnRMATCl 1HOCONSTANT s1PE16 # 4) 

WRITE (ft, l«t)XEl»XE2,xEK 

13ft FORMAT (SHITE MPE22, t /5H- NE S E?2.4/|3 H OMEGA'd,?) rEl4.4/13H 0 M EGA 
1(2*1) =E14,a/)3H 0 -BaRC 1«K) *E)4,4/t3M 0 *HaR( 2»K) =El«,4/laH O-OBA 
?R(K,2) *E13#R/13H A*BARC1 iK) =£14.4*5” (TE)/13H A*RAR( 2 ,K) =E14,4/ 
3 1 3 H P-RAR(1,2) sE14,4/13H B-BAR(2*1) »E 14.4/15H A-pRimE(?* 1 ) «E1?, 
4 4 / 1 3 H ALPHA ( K , 2) =E14,4/12H BETA(K,2) =ElS,4) 

1 a l F0RMAT(6H0N1* S1PE21.7/6H N2* »E2l,7/6H MK* *E2l,7) 

WRITE (6,146) (Z(D,Bj(n ,EA(I) ,EBCI) ,I = 1 ,NZ) 

1 4ft FORMAT (1 HO, 12X, 1 HZ, IPX, IHJ, 1 8X , 2HE A , 1 6 X , 2H£ B/ / { 1 P4 F 20 . 7) ) 

NT = N? 

NT M = VT - 1 
READ (5*t6)YA,YR 
FV = 1, 0/3,0 
YMAX = T(NT)**EY 
JV * (YMAX - 1,0) /DY 
NV = JY + 2 
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las* 

U6* 
la 7* 
la8* 
la 9* 

iso* 
151* 
152* 
1 S3* 
l sa* 

1 55* 
1?6* 
157* 
ISA* 
159* 
160* 
161 * 
162* 
163* 

1 6a* 
165* 
166 * 
167* 
16»* 
169* 
170* 
171* 

1 72* 
173* 

1 7a* 

175* 

1 76* 
177* 
17 A* 
179* 
180 * 
181* 

1 8 2 * 
183* 
lea* 
185* 
186* 
187* 

1 88 * 
189* 
190* 
191* 
192* 
193* 
i aa* 
195* 
196* 
197* 
198* 

1 99* 
200 * 
201* 
202 * 


Y(l) * 1,0 
DO 166 Js2,NY 
Y J * J - 1 
YfJ) a Y ( l ) ♦ y J*DY 
166 CONTINUE 

DO 171 Is t ♦ NY 
DO 171 Jsl,NY 
FVCI*J) a 0*0 
FTY(I*J) a 0*0 
1 7 1 CONTINUE 

DO 176 1*1, NY 
Fy(I#1) *1,0 
FIY(X»1) a t,o 
176 CONTINUE 

Do 181 J*2» NY 
JM 8 J v | 

DO 181 1*1, JM 


YY a (y(I) 
FY(ItJ) a 
F J Y ( I » j) = 
181 CONTINUE 


** 1,0) /(YfJ) * 
(1,0 * YY)*(l,o 
FV{ | f J) 


1 , 0 ) 

- Y8*YY) 


186 


191 


196 


IP * N0S0NlfPlY»X»NY,t M A)() 

IP (IE,EQ,o) GO TO 486 
GY(l) * Y (NY) . l,o 
DO 186 Js2*NY 

M(5tjwe°' 5 * (v( ' ,) * 

DO 191 Ks 1 , NY 

*Y(K) - o ,0 

DO 191 J*l,NY 

* WW>*MY(J.«> 

* 96J f y ** J *GV(I) , AY (I) ,1*1, NY) 

FORMAT (lH0,nx,t H Y f i«x, ? HGP,18)f t 2 H A Y/( lP3E 2 0. 8)) 


MO 


201 

206 

211 


216 

221 


226 


231 


J 1*1 * NT 

EX * Y(I)*Tm 

RP(I) s aY(I)*EXP(**Ey)/Y(I) 

R * R > RP ( I ) 

CONTINUE 

U R i T !, ( ^ ?06) ^ 99 ( J) » J* 1 , NY) * R 

Format (3horp//( ip E so 7)) 

*RlTE (6,211) 

FoRMAT(6M0FP(Y) ) 

DO 216 1*1, NY 

co*TiNue' J * ) 
nrITE (6,221) 

FoPMATdlHOFP INVERSE) 

Do 226 Is 1 , NY 

Do 231 1 * 1 , NY 
On 231 Jal ,NY 
EF(IfJ) * 0,0 
DO 231 Ks 1 , NY 

continue* 
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WRITE (6*31) 

pno* 


fin 236 I«1*NY 

<o5* 


wrItE (6*26) It(EF(I»J) t J»ltNY) 

206* 

236 

CONTINUE 

2f)7* 


NJK S HZ 

2o8* 


NiTKM a nTK * 1 

?0 6* 


DO 201 Mai* MY 

210* 


Y3 = Y**3 

?ll* 


00 201 l. s 1 * N T K 

212* 


TK(L*M) s T(L)/Y3 

?13* 

20t 

continue 

2 1 0* 


DO 2 06 III * NZ 

215* 


is* i 

2 1 6* 


IF (TU)-0.5) 206,246*251 

21 7* 

206 

CONTINUE 

215* 


Ism * is 

?l6* 


Go TO ?56 

2?0* 

251 

CONTINUE 

221* 


ISM * is - 1 

222* 

256 

CONTINUE 

223* 


kn = 1 

2?0* 

C 

outer loop on frequency — Index 

225* 


DO 016 Ks 1 , NY 

226* 


KK s K 

227* 


Do 225 1=1. NZ 

228* 


Tfl) = TK(I.K) 

226 * 

225 

CONTINUE 

230* 


NT* * NZ 

231* 


MTKM 5 NTK - 1 

232* 


WRITE (6*271)Y(K) 

233* 

271 

format (oh iy =E6.2) 

23«* 


write (6*276) (TKU.K) * I=1*NTK) 

235* 

27b 

Fn R MAT( 1H0* HX»2HTK//( 1PE20.7)) 

236* 


CALL WMAT(NZ*TCUT) 

237* 


00 200 I a 1 * N T K 

238* 


Do 200 Js1»nTk 

236* 


WkL m ( I » J) = *LM(I*J) 

200* 


wkL(I*J) = WL(I*J) 

201* 

200 

CONTINUE 

20?* 

205 

Continue 

203* 


DO 001 I=1*NZ 

200* 


on ooi JsltNZ 

205* 


Sw(I*J) = SH(I.J) * WKLMCI* J)*RP(K) 

206* 


WP(IfJ) * WP(I.J) + rON*RP(K)*WKL(lf J) 

2o7* 

ooi 

CONTINUE 

208* 


WRITE (6*006) 

206* 

006 

F0RMAT( 17H0WlJK(LAM8DA - 1)) 

250* 


DO Oil Is U NTK 

251* 


WRITE (6,26)I*(WKLm(T»J) *J=1*NTK) 

252* 

Oil 

continue 

253* 

016 

continue 

rv 50* 


WRITE (6,021) 

^5* 

021 

F0 R MAT(7H0S(I*J)) 

256* 


On 026 I m 1 * N T K 

257* 


WRITE (6,26)1, (SWCI.J) ,J*1, NTK) 

258* 

026 

CONTINUE 

256* 


wrITF (6*027) 

260* 

027 

FoR M AT(8H0wP(I, J) ) 



261* 


2 6?* 


263* 

428 

264* 


265* 


266* 


267* 


26ft* 


269* 

«3l 

270* 


271* 

436 

272* 


273* 

44 1 

274* 


275* 


276* 

446 

277* 


27** 


279* 


280* 


281* 


282* 

451 

283* 


284* 


285* 


286* 


287* 


288* 

456 

289* 


290* 


291* 


292* 


293* 

457 

294* 


295* 

461 

296* 


297* 


298* 


299* 


300* 

466 

301* 

471 

302* 

476 

303* 

481 

304* 


305* 

486 

306* 

491 

307* 


308* 


309* 


310* 

496 

311* 


312* 


IAGNOSTICS 


l6o 


00 42ft I a itNTK 

wrITE ( 6*26) It(WPCl» j) ♦ Jsl t NTK) 

CONTINUE 
NTKP a NIT K ♦ 1 
00 436 IbI.nTK 
00 431 Ja t . NTK 

EM(I»J) a •SW(I # J)/(r*FA(I)) 

IF(I.EO.J) E*(I»J) a FM(I.J) ♦ i.O 
CONTINUE 

EM(IfNTKP) * EBCD/EaCH 

continue 

WRITE (6.441) 
fORMAT(7HOM(I.J)) 

DO 446 1*1 .NTK 

WRITE (6,26)I.CEM(I, j) , J sl ,NTKP) 

CONTINUE 
N * NTK 
MC * NTKP 
V(t) * 4 

CALL GjRtEM.NCf NR.N.MCf 5476.JC.V) 

WRITE (6f451)(EMClfNTKP)*IalfNTK) 

FORMAT(aHOl/B//(lPE2o*7)3 
on 056 IaltNTK 
Pn * P2KI) ♦ P2* 

8(1)- l »0/EH( I f NTKP) 

ft K 1 15 * XEK/CXEl*B(in 

ft 1 C I > a (P12CI) ♦ RKfn*PK?)/PN 

CONTINUE 

DO 457 I = l * NTK 

RAU) s 0.0 

on 457 J 8 IfNtK 

RA(I) a RA(I) ♦ WP(I,J)/B(J) 

CONTINUE 

WRITE (6*461) (ZCl) f B(I) *R1(I) ♦RK(I)*RA(I)»X*1 # NTK) 

FORMAT(lHOf 13X»1HZ»|RX* 1HB* 17X»5HN?/Nl f 15X f 5HNK/Nl f ifey,3HRlK//(lP5 
E20.4)) 

PUNCH a66t CRKJ) *J*l,NTK) 

PUNCH 466* (RK( J) tJal.NTK) 

PUNCH 466* (RA(J) *J*1 ,NTK) 

FnRMATC5E15.6) 

call exit 
write ( 6 f 48 1 ) 

FOBMATCSIHOERROR IN go, SOLVING) 

Go TO 471 

WRITE (6,491)IE*NY 

FORMAT ( ISHOERROR in FP(Y)/5H if al3.5X.4HNy *13) 

WRITE (6*211) 

DO 496 la 1 » NT 

WRITE (6*26)1* (FY( I. J) * Jsl » NT) 

CONTINUE 
Go TO 471 
End 


5.63 CPU SECONDS 


ION T IMF * 



l6l 


1* 

?* C 

3* 

a* 

s* 

6 * 

7* 

B* 

9* 

1 0* 

1 1* 

\ ?* 
ti* 

ta* 

15* 
l a* 

17* 

1 A* 

1** 

?0* 


SH&RO"Tlht WMAT(NZ#TCU7) ’ 

CALCULATES W(lAMBDA) AND W(L4M0DA - 1) at L?NE center 
DOUBLE PRECISION T»El *t2fE3.E4 t Al*A2tA3»A«» A5.ESB.pt » r?»r3 

double precision f» ft . xk 

C0MH0N/BLK1/T(31) ♦El(3l)»E2(31)»Ei(3l) .E«(3t) 

COMMOW/Bl K2/F(3tf 31) f FIC3l*31) »WL(31f 31) fWLMCJl, 31) tXKC3l ♦3t)« 
IX L (31»31) 

10 PO«MAT(IH0*I4,2X»1P5e?O.7/(7X«*;E2O,7)) 

NT * N? 

DO t5 |*liNZ 
DO 15 J*lfNZ 
*L(I*J) * 0.0 
i*L m ( I ♦ J) » 0.0 
15 CONTINUE 
NM a NZ m 1 
Do Js 1 » NM 

JP s I + 1 
DO 20 JSJP.NM 
A 1 = DABS(T( ) 

IF fAl.G7.TCUT) &G Tn 25 



162 


21* 
22 * 
23* 
26* 
25* 
26* 
27* 
20 * 
2^* 
30* 
3 l * 
32* 
33* 
3«* 
3*5* 
36* 
37* 
36* 
39* 

ao* 

ai* 

a2* 

63 * 

66* 

6*5* 

66 * 

67* 

40* 

46* 

50* 

51* 

52* 

53* 

56* 

55* 

56* 

57* 

58* 

59* 

60* 

61* 

62* 

63* 

66 * 

65* 

66 * 

67* 

68 * 

69* 

70 * 

71* 

72* 

73 * 

76* 

75* 

76* 

77* 

78* 


A2 
A3 
A6 
A5 
81 
1 ♦ 
B2 
1 - 


s 0A8SCT(J) - TCI)) 

* OaBS(T(J+J) - TfD) 
a TfJ) - T(J-1) 

* T ( J+ 1 ) - T(J) 

* (T( I)*(ESB( Al - ES0CA2»2)) 
tSB{A2*3) • ES9(Aj f 3))/A6 

* (T(J»1)/A6)*(ESbCA 1*2) - FSB(A2*2)) - 
E$B( A3f 2) ) 


♦ DEXP(-Al) » OrXPC- 


C T ( J* 1 ) /A5) * 


93 * CT(I)*CtS0(Aa f 2) . ESB( A3» 2) ) 
1 ♦ ES8( A3*3) - ESb( A? ,3))/a5 
* O.5*(01 „ m m bj) 
WLHCIiJ) * WL(I,J) 

20 CONTINUE 
25 CONTINUE 


♦ OEXPUA?) 


OEXP(« 


A?) 

CESb(A2»?) 

*3) 


00 35 J*2 * NM 
IP « J ♦ 1 


DO 30 I«IP»NZ 


A3 » DaBS(TCI) * T( J+ l ) ) 

IF (A3,GT*TCUT) Go To 35 
Al S OA0S(T(I) - T(J-D) 

A 2 * DABSCT(I) - T ( J) ) 

A6 ■ T ( J) - T(J*1) 

AS « T( J* 1 ) - T(J) 

Bl » (Al/A6)*(E$B(A2«2) ^ ESB( 
1 ESB(A3.2) * ESB( A2» 2) ) 

B2 a (0EXP(«A3) » DExP("*A2) ♦ 
83 a (OEXP (*• A2) • DEXP(-Al) ♦ 
WL(I»J) a 0,5*(B1 ♦ R2 - R3) 
WlNCItJ) * Wl(IfJ) 


Alf2)) ♦ C(T(J+1) . t(I))/A5>*( 

ES6(A?*3) - tSB(A3*3i)/A5 
ESB ( A 1 *3) - ESB(A?,3) )/A6 


30 CONTINUE 
35 CONTINUE 

00 60 Ia2»wM 

Al a T(I) * t(MJ 

A2 a TCl+t) • TCI) 

01 a 2,0 - ESB( Al *2) - ESB(A2*2) 

B2 « (0.5 * 0EXP(*A1) * ESB ( A 1 * 3) ) /A 1 

03 * (0,5 * 0EXP(»A2) ♦ ES8( A2 t 3) ) /A? 

*L ( I ♦ I) a 0,S*(B1 m R2 * 83) 

WLH(I,I) a WL(I,I) - 1,0 
60 CONTINUE 


DO 65 lai,NM 

Al a DAB5(T(NM)-T(D) 

IF (Al,GT f TCUT) GO To 65 
A2 a DAB$(T(NZ) * T C T) ) 

A3 a T(NZ) - T(NM) 

92 * (T(NM)/A3)*(ESB(A1»?) * EsS(A 2 ,?)) 

WL(IiNZ) i Bl , B2 
«lH(i*NZ) a «UI»NZ) 

65 CONTINUE 

A3 a T(2) - TCI) 

DO 50 Ia2*NZ 

A2 a DABS(T(2) * TCI)) 

IF (A 2 ,GTtTCUT) Go j 0 50 
Al a DABS(T(1) - T(I)) 

81 * ((T(2) - T(I))/a3)*(ESB(A2»2) - ESB(A1,?>) 

52 a (0FXP(-A2) - OExP(-Al) ♦ ESB(A1,3) . ESP ( A2* 3) ) /A3 



163 


79* 

|*0* 

M* 

a?* 

*3* 

ao* 

«5* 

a6* 

a7* 

aa* 

rp* 

90* 

91 * 

9?* 

93* 

qa* 

95* 

96* 

97* 

qfi* 

99* 

100* 

101 * 

102* 

10 3* 

104* 

105* 

agnostics 


w L ( I * 1 ) c 0,5* (81 + R2) 

WL.M{I»1) * HUM) 

50 CONTINUE 

TF (A3.GT.TCUT) GO Tn 55 

w l { 1 1 1 ) a 0.5*(1,0 - E$8(A3*2) «- Co.S - OtxP(-A3) + 
GO TO 60 

55 wtUfll = 0,5 - (0,25 /a3) 

60 |) a WLCl'tO - 1.0 

A 1 s T(NZ) - T(N^) 

IF -(■•Al-.GT.TCUT) &Q To 65 

wtCNZ»MZ3 a 0,5* (1,0 * ES8(Al»?) - (0,5 - DEXP(-Al) 
1 /Al) 

GO TO 70 

65 HCNZtNZ) * 0,5 - (0.25/A1) 

70 ^lN(NZ,NZ) = WL,(NZ,N7) * 1,0 

- wrITE (6.75) 

75 FORMAT ( 10H0W (LAH0DA) ) 

DO 80 I * t » N T 

WRITE (6,10) I,(Wt(I»J) ,JalfNT) 

ao Continue - 

WRITE (6,85) 

85 F0«MaT( UHOWCLaNBDA - 1)) 

Do 90 lei f NT 

WRITE (6,10) I,(WLH(I # J) ,ja l# NT) 

90 CONTINUE 
RETURN 
&Nf> 


E$8(A3,3))/A3) 


+ E S P ( A 1 * 3 ) ) 


ON TIME a 


2,60 CPU SECONDS 
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1 * 

2 * 

3 * 

a* 

5 * 

6 * 

7 * 

8 * 

9 * 

[ 0 * 

1 i* 
12 * 
13 * 
19 * 

1 5 * 

16 * 
j 7 * 
t 8 * 
19 * 
20 * 
21 * 
22 * 
? 3 * 

2 «* 

25 * 


r 

r 

r 

c 

c 

c 

c 

c 

c 


SUBROUTINE GJR(A*NC»NR» N»MC»J *JC» V) 
DOUBLE PRECISION A»X,V 
dimension acnr,nC),jc(1)»v(2) 


jc is the permutation vector 

KO IS the OPTION key for DETERMINANT EVALUATION 
KI Is THE OPTION KEY FOR MATRIX INVERSION 
L IS THE COUJHN CONTROL f OR AX*B 

m is the emu** contol for matrix inversion 


INITIALIZATION 


iweV(l) 

M* 1 

SsU 

LpN^(MC*N)»( IN/R) 
KD«2-MOD( iw/2f 2) 
IF(KO.EQ.t) V(?)aO, 
KI#2-MO0{ IW»2) 

GO TO { 1 0*20) *KI 


initialize JC for inversion 


C 

C 

C 


10 DO 15 I* 1 1 N 

15 jcn)*i 
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|P6* 

It* 

;>8* 

?4* 
"SO* 
31* 
3?* 
33* 
34* 
35* 
36* 
3 7* 
36* 
3^* 
40* 
41 * 
4?* 
43* 

4 4* 

45* 

46* 

47 * 
46 * 
44 * 
50 * 
51* 
52* 
53* 
54* 
55* 
56* 
57* 
58* 
54 * 
* 0 * 
61 * 
62 * 
6 3* 
64* 
65* 
66 * 
67* 
68 * 

6 ** 

70* 

7!* 

72* 

73* 


74* 
75* 
76* 
V7 7* 

#7«* 

74* 

80 * 

61* 


62 * 

A3* 


c — - 

r sf arch for pivot row 

c — — 

2 0 on 86 ?* t » N 

GO TO ( 30 • 25 ) » k I 
25 Mel 

30 IF (I,FQ*M) GO TO 55 

Xs* 1 • 

DO 35 J=I»N 

IF (X.gT,ABS(A( J*IJ) j GO TO 35 
X a DABS ( A ( J * 1 ) } 

KsJ 

35 Continue 

if (k.fq.D go to 55 

s*-s 

V(l)s*V(t ) 

GO TO (40*45) *KI 
40 Mysjcd) 

JC(I)*JC(K) 

JC(K)«MU 

C 

c interchange row i and row k 

c 

45 on 50 JSM.L 
X*A( I* J) 

A( It J) sA(K# J) 

50 A(K*J)sX 

C 

r test for singularity 

c 

5S If (OARS(A(I»I))*GT.O.O) GO TO 60 

C — 

C MATRIX IS SINGULAR 

r. * — «* — --- 

lF(KD,FQ,n V( 1 )« 0 . 

JC(t)*I*l 

RETURN 6 

60 GO TO (65»70) t KD 

C 

c compute the determinant 

c — •* — 

65 3F(A<I.I).LTtO.) S»«*s 

V(2) s V (2) + OLOG(DABS(A(I»I))) 

70 Xs A ( I * J ) 

Afl* I)al, 

C REDUCTION Of THE I*TM RQW 

C — 

on 75 JsM» L 
A(I.J)s4(I,J)/X 

C 

c TEST OVERFLOW SWITCH, IF ON 

0 return negative VALUE OF I IN <JC ( 1 ) 

0 ... »*• 

CALL dverfl (ifl) 

IF (IFL.EO.l) GO TO 120 
75 Continue 

o 
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04* 
05* 
06* 
0 7* 
6 0 * 
09 * 
90* 
9t* 
9 ?* 
93 * 
94* 
95* 
96* 
97* 
90* 
99* 
10 0 * 
101 * 
to?* 

10 3* 
1 0 4* 
105* 
106 * 
1 07* 
100 * 
109* 
110 * 
111 * 
11 ?* 
1 1 3* 
114 * 
115* 
116* 
1 t 7* 

1 1 8 * 
119* 
5 ? 0 * 
1 21 * 
1 ??* 
1?3* 
1?4* 
1?5* 


c REDUCTION Of all REMAINING rows 

C * 

Of) 65 K*1*N 
IF (K.FG.I) GO TO 65 
X a A ( K * I ) 

A(K»I)*0# 
on 60 J*MfL 

A(K, J)«A(K.J)-X*A(I,j) 

C - 

r TEST OVERFLOW SWITCH. IF ON 

C RETURN NEGATIVE VALUE Of I IN JC(1) 


call overfl (ifl) 

IF ( IFL.FO. n GO TO 1?0 
00 CONTINUE 
85 CONTINUE 

c — 

f AX80 AND DEt.(A) are NOW COMPUTED 


GO TO (9o*U5) ,KI 

C — 

c PERMUTATION of the COLUMNS FOR matrix INVERSION 


90 DO no Jsl.N 

IF (JCCJ).EQ*J) GO TO 110 

DO 95 I»JJ»N 

IF CJC(I).E0#J) GO To 100 
95 CONTINUE 

too jccn*jc(j) 

DO 105 Kal.N 
X«A(K*I) 

A(K*I)aA(K,J) 

105 A ( K » J) aX 
110 CONTINUE 
115 JCCIJsN 

IF(KD.EO.I) V ( 1 ) *S 

return 

1?0 JC(l)«l"I 

IFCKD.EO.l) V(i)*s 

RETURN 6 

END 


I AGNOSTICS 

TION TIME * 1.96 CPU SECONDS 
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1* 

2* 

3* 

a* 

5* 
6 * 
7* 
8 * 
9* 
1 0* 
1 1* 
12* 
1.3* 
1 4* 
IS* 
16* 
1 7* 
lfi* 
1** 
20 * 
21 * 
22 * 
23* 

24* 
25* 
26 * 
27* 
? 8 * 
?9* 
30* 
31* 
32* 
33* 
3 0* 


nnUAU PRECISION FUNCTION NOSONIC A.x.L.LMAk) 
DOUBLE PRECISION A»X,F 
DIMENSION A(l).XCl) 

M * L - 1 
MAX * Ni*LMAX ♦ L 
MAX s N*L*AX + L 
DO 10 1*1. 1 
X(I) 8 1.0 
10 CONTINUE 
K 1 * * LM AX 
00 55 K*1 »l 
K 1 * K 1 ♦ UMAX 

IF *A(k2) ) *15.80. 15 
IS Do 30 T * 1 » L 
Jl ■ K1 ♦ I 
IF ( AC Jl ) ) 20 ♦ 30 » 20 
20 F * I.O/ACJD 

xcn * x c i ) * f 

00 ?$ JlsI.MAX.LMAX 
A f J 1 ) « A ( J 1 ) * F 
25 CONTINUE 
30 CONTINUE 

A ( K 2 ) s X(K) 

X(K) = 1.0 
DO 50 1*1. L 
KJ « K - I 
IF (KI) 35,50.35 
35 Jl ■ Kt ♦ I 

IF ( A C J 1 ) ) 40.50.40 
40 A C J 1 ) s 0.0 

DO a5 J2sI.MAX.tMAX 
Jl = J2 + Kl 
A ( J?) = A ( J2) - A(Jn 
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35* 

as 

CONTINUE 

56* 

so 

Continue 

57* 

ss 

continue 

38* 


00 70 !»1»M 

39* 


IF (X(D) 60»80,60 

#0* 

60 

F * 1,0/XCI) 

at* 


DO 65 JlsI.MAXtlMAX 

a?* 


A(Jt) s A(J1)*F 

a3* 

65 

CONTINUE 

aa* 

70 

continue 

as* 


nosoni » 1 

a6* 

75 

return 

a7* 

00 

NOSONI B 0 

ae* 


GO TO 75 

a9* 


END 

AGNOSTICS 

ION TIMS 

X 

,98 CPU SECONDS 
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1* 


DOUBLE PRECISION FUNCTION ESB(X*N) 

2* C 


CALCULATES E2 * E3» E4. FROM El 

3* 


double precision x 

a* 


DOUBLE PRECISION EXlN 

5* 


IF (X.LE.0.0) GO TO 35 

ft* 


E * EXIN(X) 

7* 


GO TO ( 25* 10* 10* 10) *K) 

a* 

10 

EX s OEXP(-X) 

9* 


E a EX * X*E 

10* 


Go TO f 25»?5» 15t 15) *N 

ll* 

15 

F a 0 .5* ( EX “ X*E) 

1?* 


GO TO (25»25*25t*0 ) *m 

M* 

20 

F = (Ex - X*E)/3.0 

14* 

25 

ESB = F 

i 5* 

30 

RETURN 

1 ft* 

35 

GO TO (40*45*50*55) iM 

1 7* 

40 

ESB = 0.0 

i a* 


Go TO 30 

1 9* 

45 

ESB « 1 *0 

20* 


GO TO 30 

?1* 

50 

EsB * 0.5 

22* 


GO TO 30 

23* 

55 

ESB « l. 0/3.0 

24/* 


GO TO 30 

?S* 


end 

IAGNOSTICS 

TfpN T I Mp 

s 

,4ft CPU SECONDS 
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1* 

?.* 

3* 

a* 

5* 

6 * 

7* 

8 * 

9* 

10 * 

n* 
1 2 * 
J3* 

t** 

1*5* 

16* 

17* 

1** 

19* 

? 0 * 

?i* 

22 * 

23* 


DOUBLE PRECISION FUNCTION EXIN(Y) 

CALCULATE S El 
DOUBLE PRECISION Y 

double precision x» a*b?c»e*bb»cc 

DIMENSION A{6) , B<4) ,C(4) 

* 9 • 9999 1^30- 01 »• 2 .4991055D" 01 *5.51 99680*0? • 
l»9«76004D+03f 1*078570*03/ 

DATA B/8,5733287^01D+O0f 1.8059 0 1697300+0 l»8,fe34760892*5D+00 t 

12.6777373430+01/ 

DATA C/9, 57332234540+00 .2. 5632956 14860+01. 2, i 09965x0*270+0 1. 
13.95849692280+00/ 

X » Y 

IF (X+1 ,0) 10+15+15 

E * All) + X*(A(2) + X*( A(3) + X* ( A ( 4) + X*CA{5) ♦ X*A(6)))>) 

E * E m DLOG(X) 

EXIN » E 
RETURN 

15 Bb B 8(4) + X*ce(3) + X* ( B ( 2 ) + X*(B(t) ♦ X))) 

CC » C(4) ♦ X*(C(3) + X*(C(2) + X*(C(l) + X))) 

E • (BB/CC)*DEXP(+X)/X 

EXlN ■ E 

RETURN 

end 


10 


AGNOSTICS 
ION TIME 


,46 CPU SECONDS 


TPFS M 80NF2J, 
024+08/30*23104 
6 REL 
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1* DOUBLE PRECISION FUNCTION ESBCX.N) 

2* C CALCULATES E2» ES* E4» FROM El 

3 * Double precision x 

a* double precision exin 

5* IP (X.LE.0.05 00 to 35 

6* E P EXIN(X) 

7* 60 TO (25»10M0»I0)»n 

8* 10 Ex p D|XP<*X) 

?f e i EX • X*E 

lor 00 TO <25t25M5.1S5,N 

iir is e p o.ittix • x*i> 

lif 60 TO (25f 25*25.20) «N 

Iff 20 E P (EX * XfC)/3,0 

l«f 25 EBB f | 

l5f 30 RETURN 

i** 55 60 TO («0*45.50*SS).ki 

I?f ao til * o.o 

is* 60 TO 30 

Iff «5 EBB P 1*0 

20* 60 TO 30 

21* 30 EBB P 0.5 

gl* 60 TO 30 

23f 55 EsB p I. 0/3.0 

gap 60 TO 30 

25f end 


oi agnostic® 

1TI0N TIME • ,«1 "CPU SECONDS 



double PREClflON FUNCTION EXINCY) 

CALCULATES El 

double precision y 

DOUBLE PRECISION X, A.BfC.E.BB.CC 
DIMENSION A(S)fB(«).C(A) 

DATA A/»5'.772l56D*01, 9. 99991930*01 ♦•2*A991055 Dp 01 »5.5!9968D*02i 
1*9.?60 P'Dp0T* I ; 0TO37D.0I/ 

DATA B/S .5T332B7S01DP 00* i .605*0 16*750 D+01 .8.63 07608925D+00. 

12* 6777373630*01/ 

DATA C/9" 5733223656D+00* 2, 5632936 1666DT01 .2,1 0996530827D701 * 

13. 95649692260*00/ 

X p Y 

ip tx*t,-o r roiftrs : — ; — 

10 E » AC II ♦ X*( A(2) ♦ x* ( A ( 35 ♦ X»(A(A) ♦ X*(A(5) 4 X*A(65)5)J 
E P E » DL06(X) 

EXIN p C 
RETURN 

15 BB P BCJ) ♦ X» (B(3 ) ♦ X*(B{2) 4 X*(B(t) 4 X))) 

CC p 1T(«T ♦ >T*(C(3) + X*(C(2) 4 X*(C( 1) 4 X5 ) > 

E p (BB/CC)*DEXP(*X)/X 
EXIN p E 

return 

end 


I AGNOSTICS 


TION TIME * 


,64 CPU SECONDS 



